text:
module Plotting
import Luna: Grid, Maths, PhysData, Processing
import Luna.PhysData: wlfreq, c, ε_0
import Luna.Output: AbstractOutput
import Luna.Processing: makegrid, getIω, getEω, getEt, nearest_z
import PyPlot: ColorMap, plt, pygui, Figure
import FFTW
import Printf: @sprintf
import Base: display

"""
    displayall()

`display` all currently open PyPlot figures.
"""
function displayall()
    for fign in plt.get_fignums()
        fig = plt.figure(fign)
        display(fig)
    end

end

display(figs::AbstractArray{Figure, N}) where N = [display(fig) for fig in figs]

"""
    cmap_white(cmap, N=512, n=8)

Replace the lowest colour stop of `cmap` (after splitting into `n` stops) with white and
create a new colourmap with `N` stops.
"""
function cmap_white(cmap; N=2^12, n=8)
    vals = collect(range(0, 1, length=n))
    vals_i = collect(range(0, 1, length=N))
    cm = ColorMap(cmap)
    clist = cm(vals)
    clist[1, :] = [1, 1, 1, 1]
    clist_i = Array{Float64}(undef, (N, 4))
    for ii in 1:4
        clist_i[:, ii] .= Maths.BSpline(vals, clist[:, ii]).(vals_i)
    end
    ColorMap(clist_i)
end

"""
    cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)

Make an array of `num` different colours that follow the colourmap `cmap` between the values
`cmin` and `cmax`.
"""
function cmap_colours(num, cmap="viridis"; cmin=0, cmax=0.8)
    cm = ColorMap(cmap)
    n = collect(range(cmin, cmax; length=num))
    cm.(n)
end

"""
    subplotgrid(N, portrait=true, kwargs...)

Create a figure with `N` subplots laid out in a grid that is as close to square as possible.
If `portrait` is `true`, try to lay out the grid in portrait orientation (taller than wide),
otherwise landscape (wider than tall).
"""
function subplotgrid(N, portrait=true; colw=4, rowh=2.5, title=nothing)
    cols = ceil(Int, sqrt(N))
    rows = ceil(Int, N/cols)
    portrait && ((rows, cols) = (cols, rows))
    fig, axs = plt.subplots(rows, cols, num=title)
    ndims(axs) > 1 && (axs = permutedims(axs, (2, 1)))
    if cols*rows > N
        for axi in axs[N+1:end]
            axi.remove()
        end
    end
    fig.set_size_inches(cols*colw, rows*rowh)
    fig, N > 1 ? axs : [axs]
end

"""
    get_modes(output)

Determine whether `output` contains a multimode simulation, and if so, return the names
of the modes.
"""
function get_modes(output)
    t = output["simulation_type"]["transform"]
    !startswith(t, "TransModal") && return false, nothing
    lines = split(t, "\n")
    modeline = findfirst(li -> startswith(li, "  modes:"), lines)
    endline = findnext(li -> !startswith(li, " "^4), lines, modeline+1)
    mlines = lines[modeline+1 : endline-1]
    labels = [match(r"{([^,]*),", li).captures[1] for li in mlines]
    angles = parse.(Float64, [match(r"ϕ=([0-9]+.[0-9]+)π", li).captures[1] for li in mlines])
    if !all(angles .== 0)
        for i in eachindex(labels)
            if startswith(labels[i], "HE")
                if angles[i] == 0
                    θs = "y"
                elseif angles[i] == 0.5
                    θs = "x"
                else
                    θs = "$(angles[i])π"
                end
                labels[i] *= " ($θs)"
            end
        end
    end
    return true, labels
end

"""
    stats(output; kwargs...)

Plot all statistics available in `output`. Additional `kwargs` are passed onto `plt.plot()`
"""
function stats(output; kwargs...)
    stats = output["stats"]

    pstats = [] # pulse statistics
    haskey(stats, "energy") && push!(pstats, (1e6*stats["energy"], "Energy (μJ)"))
    for (k, v) in pairs(stats)
        startswith(k, "energy_") || continue
        str = "Energy "*replace(k[8:end], "_" => " ")*" (μJ)"
        push!(pstats, (1e6*stats[k], str))
    end
    for (k, v) in pairs(stats)
        startswith(k, "peakpower_") || continue
        Pfac, unit = power_unit(stats[k])
        str = "Peak power "*replace(k[11:end], "_" => " ")*" ($unit)"
        push!(pstats, (Pfac*stats[k], str))
    end
    if haskey(stats, "peakpower")
        Pfac, unit = power_unit(stats["peakpower"])
        push!(pstats, (Pfac*stats["peakpower"], "Peak power ($unit)"))
    end
    haskey(stats, "peakintensity") && push!(
        pstats, (1e-16*stats["peakintensity"], "Peak Intensity (TW/cm\$^2\$)"))
    haskey(stats, "fwhm_t_min") && push!(pstats, (1e15*stats["fwhm_t_min"], "min FWHM (fs)"))
    haskey(stats, "fwhm_t_max") && push!(pstats, (1e15*stats["fwhm_t_max"], "max FWHM (fs)"))
    haskey(stats, "fwhm_r") && push!(pstats, (1e6*stats["fwhm_r"], "Radial FWHM (μm)"))
    haskey(stats, "ω0") && push!(pstats, (1e9*wlfreq.(stats["ω0"]), "Central wavelength (nm)"))

    fstats = [] # fibre/waveguide/propagation statistics
    if haskey(stats, "electrondensity")
        push!(fstats, (1e-6*stats["electrondensity"], "Electron density (cm\$^{-3}\$)"))
        if haskey(stats, "density")
            push!(fstats,
                 (100*stats["electrondensity"]./stats["density"], "Ionisation fraction (%)"))
        end
    end
    haskey(stats, "density") && push!(
        fstats, (1e-6*stats["density"], "Density (cm\$^{-3}\$)"))
    haskey(stats, "pressure") && push!(
        fstats, (stats["pressure"], "Pressure (bar)"))
    haskey(stats, "dz") && push!(fstats, (1e6*stats["dz"], "Stepsize (μm)"))
    haskey(stats, "core_radius") && push!(fstats, (1e6*stats["core_radius"], "Core radius (μm)"))
    haskey(stats, "zdw") && push!(fstats, (1e9*stats["zdw"], "ZDW (nm)"))

    z = stats["z"]*1e2

    multimode, modes = get_modes(output)

    Npl = length(pstats)
    if Npl > 0
        pfig, axs = subplotgrid(Npl, title="Pulse stats")
        for n in 1:Npl
            ax = axs[n]
            data, label = pstats[n]
            multimode && (ndims(data) > 1) && (data = data')
            ax.plot(z, data; kwargs...)
            ax.set_xlabel("Distance (cm)")
            ax.set_ylabel(label)
            multimode && (ndims(data) > 1) && ax.semilogy()
            multimode && (ndims(data) > 1) && ax.legend(modes, frameon=false)
        end
        pfig.tight_layout()
    end

    Npl = length(fstats)
    if Npl > 0
        ffig, axs = subplotgrid(Npl, title="Other stats")
        for n in 1:Npl
            ax = axs[n]
            data, label = fstats[n]
            multimode && (ndims(data) > 1) && (data = data')
            ax.plot(z, data; kwargs...)
            ax.set_xlabel("Distance (cm)")
            ax.set_ylabel(label)
            multimode && (ndims(data) > 1) && should_log10(data) && ax.semilogy()
            multimode && (ndims(data) > 1) && ax.legend(modes, frameon=false)
        end
        ffig.tight_layout()
    end
    [pfig, ffig]
end

"""
    should_log10(A, tolfac=10)

For multi-line plots, determine whether data for different lines contained in `A` spans
a sufficiently large range that a logarithmic scale should be used. By default, this is the
case when there is any point where the lines are different by more than a factor of 10.
"""
function should_log10(A, tolfac=10)
    mi = minimum(A; dims=2)
    ma = maximum(A; dims=2)
    any(ma./mi .> 10)
end

window_str(::Nothing) = ""
window_str(win::NTuple{4, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win[2:3]...)
window_str(win::NTuple{2, Number}) = @sprintf("%.1f nm to %.1f nm", 1e9.*win...)
window_str(window) = "custom bandpass"

"""
    prop_2D(output, specaxis=:f)

Make false-colour propagation plots for `output`, using spectral x-axis `specaxis` (see
[`getIω`](@ref)). For multimode simulations, create one figure for each mode plus one for
the sum of all modes.

# Keyword arguments
- `λrange::Tuple(Float64, Float64)` : x-axis limits for spectral plot (wavelength in metres)
- `trange::Tuple(Float64, Float64)` : x-axis limits for time-domain plot (time in seconds)
- `dBmin::Float64` : lower colour-scale limit for logarithmic spectral plot
- `resolution::Real` smooth the spectral energy density as defined by [`getIω`](@ref).
"""
function prop_2D(output, specaxis=:f;
                 trange=(-50e-15, 50e-15), bandpass=nothing,
                 λrange=(150e-9, 2000e-9), dBmin=-60,
                 resolution=nothing, modes=nothing,
                 kwargs...)
    z = output["z"]*1e2
    if specaxis == :λ
            specx, Iω = getIω(output, specaxis, specrange=λrange, resolution=resolution)
    else
            specx, Iω = getIω(output, specaxis, resolution=resolution)
    end

    t, Et = getEt(output, trange=trange, bandpass=bandpass)
    It = abs2.(Et)

    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
    specx .*= specxfac

    multimode, modelabels = get_modes(output)

    if multimode
        fig = _prop2D_mm(modelabels, modeidcs(modes, modelabels), t, z, specx, It, Iω,
                         speclabel, speclims, trange, dBmin, window_str(bandpass);
                         kwargs...)
    else
        fig = _prop2D_sm(t, z, specx, It, Iω,
                         speclabel, speclims, trange, dBmin, window_str(bandpass);
                         kwargs...)
    end
    fig
end

modeidcs(m::Int, ml) = [m]
modeidcs(m::Symbol, ml) = (m == :sum) ? [] : error("modes must be :sum, a single integer, or iterable")
modeidcs(m::Nothing, ml) = 1:length(ml)
modeidcs(m, ml) = m

# Helper function to convert λrange to the correct numbers depending on specaxis
function getspeclims(λrange, specaxis)
    if specaxis == :f
        specxfac = 1e-15
        speclims = (specxfac*c/maximum(λrange), specxfac*c/minimum(λrange))
        speclabel = "Frequency (PHz)"
    elseif specaxis == :ω
        specxfac = 1e-15
        speclims = (specxfac*wlfreq(maximum(λrange)), specxfac*wlfreq(minimum(λrange)))
        speclabel = "Angular frequency (rad/fs)"
    elseif specaxis == :λ
        specxfac = 1e9
        speclims = λrange .* specxfac
        speclabel = "Wavelength (nm)"
    else
        error("Unknown specaxis $specaxis")
    end
    return speclims, speclabel, specxfac
end

# single-mode 2D propagation plots
function _prop2D_sm(t, z, specx, It, Iω, speclabel, speclims, trange, dBmin, bpstr; kwargs...)
    num = "Propagation" * ((length(bpstr) > 0) ? ", $bpstr" : "")
    pfig, axs = plt.subplots(1, 2, num=num)
    pfig.set_size_inches(12, 4)
    Iω = Maths.normbymax(Iω)
    _spec2D_log(axs[1], specx, z, Iω, dBmin, speclabel, speclims; kwargs...)

    _time2D(axs[2], t, z, It, trange; kwargs...)
    pfig.tight_layout()
    return pfig
end

# multi-mode 2D propagation plots
function _prop2D_mm(modelabels, modes, t, z, specx, It, Iω,
                    speclabel, speclims, trange, dBmin, bpstr;
                    kwargs...)
    pfigs = Figure[]
    Iω = Maths.normbymax(Iω)
    id = "($(string(hash(gensym()); base=16)[1:4])) "
    for mi in modes
        num = id * "Propagation ($(modelabels[mi]))" * ((length(bpstr) > 0) ? ", $bpstr" : "")
        pfig, axs = plt.subplots(1, 2, num=num)
        pfig.set_size_inches(12, 4)
        _spec2D_log(axs[1], specx, z, Iω[:, mi, :], dBmin, speclabel, speclims; kwargs...)

        _time2D(axs[2], t, z, It[:, mi, :], trange; kwargs...)
        push!(pfigs, pfig)
    end

    num = id * "Propagation (all modes)" * ((length(bpstr) > 0) ? ", $bpstr" : "")
    pfig, axs = plt.subplots(1, 2, num=num)
    pfig.set_size_inches(12, 4)
    Iωall = dropdims(sum(Iω, dims=2), dims=2)
    _spec2D_log(axs[1], specx, z, Iωall, dBmin, speclabel, speclims; kwargs...)

    Itall = dropdims(sum(It, dims=2), dims=2)
    _time2D(axs[2], t, z, Itall, trange; kwargs...)
    pfig.tight_layout()
    push!(pfigs, pfig)

    return pfigs
end

# a single logarithmic colour-scale spectral domain plot
function _spec2D_log(ax, specx, z, I, dBmin, speclabel, speclims; kwargs...)
    im = ax.pcolormesh(specx, z, 10*log10.(transpose(I)); shading="auto", kwargs...)
    im.set_clim(dBmin, 0)
    cb = plt.colorbar(im, ax=ax)
    cb.set_label("SED (dB)")
    ax.set_ylabel("Distance (cm)")
    ax.set_xlabel(speclabel)
    ax.set_xlim(speclims...)
end

# a single time-domain propagation plot
function _time2D(ax, t, z, I, trange; kwargs...)
    Pfac, unit = power_unit(I)
    im = ax.pcolormesh(t*1e15, z, Pfac*transpose(I); shading="auto", kwargs...)
    cb = plt.colorbar(im, ax=ax)
    cb.set_label("Power ($unit)")
    ax.set_xlim(trange.*1e15)
    ax.set_xlabel("Time (fs)")
    ax.set_ylabel("Distance (cm)")
end

"""
    time_1D(output, zslice, y=:Pt, kwargs...)

Create lineplots of time-domain slice(s) of the propagation.

The keyword argument `y` determines
what is plotted: `:Pt` (power, default), `:Esq` (squared electric field) or `:Et` (electric field).

The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.

The keyword argument `oversampling` determines the amount of oversampling done before plotting.

Other `kwargs` are passed onto `plt.plot`.
"""
function time_1D(output, zslice=maximum(output["z"]);
                y=:Pt, modes=nothing,
                oversampling=4, trange=(-50e-15, 50e-15), bandpass=nothing,
                FTL=false, propagate=nothing,
                kwargs...)
    t, Et, zactual = getEt(output, zslice,
                           trange=trange, oversampling=oversampling, bandpass=bandpass,
                           FTL=FTL, propagate=propagate)
    if y == :Pt
        yt = abs2.(Et)
    elseif y == :Et
        yt = real(Et)
    elseif y == :Esq
        yt = real(Et).^2
    else
        error("unknown time plot variable $y")
    end
    multimode, modestrs = get_modes(output)
    if multimode
        if modes == :sum
            y == :Pt || error("Modal sum can only be plotted for power!")
            yt = dropdims(sum(yt, dims=2), dims=2)
            modestrs = join(modestrs, "+")
            nmodes = 1
        else
            isnothing(modes) && (modes = 1:length(modestrs))
            yt = yt[:, modes, :]
            modestrs = modestrs[modes]
            nmodes = length(modes)
        end
    end

    yfac, unit = power_unit(abs2.(Et), y)

    sfig = plt.figure()
    if multimode && nmodes > 1
        _plot_slice_mm(plt.gca(), t*1e15, yfac*yt, zactual, modestrs; kwargs...)
    else
        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
        label = multimode ? zs.*" ($modestrs)" : zs
        for iz in eachindex(zactual)
            plt.plot(t*1e15, yfac*yt[:, iz]; label=label[iz], kwargs...)
        end
    end
    plt.legend(frameon=false)
    add_fwhm_legends(plt.gca(), "fs")
    plt.xlabel("Time (fs)")
    plt.xlim(1e15.*trange)
    ylab = y == :Et ?  "Field ($unit)" : "Power ($unit)"
    plt.ylabel(ylab)
    y == :Et || plt.ylim(ymin=0)
    sfig.set_size_inches(8.5, 5)
    sfig.tight_layout()
    sfig
end

# Automatically find power unit depending on scale of electric field.
function power_unit(Pt, y=:Pt)
    units = ["kW", "MW", "GW", "TW", "PW"]
    Pmax = maximum(Pt)
    oom = clamp(floor(Int, log10(Pmax)/3), 1, 5) # maximum unit is PW
    powerfac = 1/10^(oom*3)
    if y == :Et
        sqrt(powerfac), "$(units[oom])\$^{1/2}\$"
    else
        return powerfac, units[oom]
    end
end

"""
    spec_1D(output, zslice, specaxis=:λ, log10=true, log10min=1e-6)

Create lineplots of spectral-domain slices of the propagation.

The x-axis is determined by `specaxis` (see [`getIω`](@ref)).

If `log10` is true, plot on a logarithmic scale, with a y-axis range of `log10min`.

The keyword argument `modes` selects which modes (if present) are to be plotted, and can be
a single index, a `range` or `:sum`. In the latter case, the sum of modes is plotted.

Other `kwargs` are passed onto `plt.plot`.
"""
function spec_1D(output, zslice=maximum(output["z"]), specaxis=:λ;
                 modes=nothing, λrange=(150e-9, 1200e-9),
                 log10=true, log10min=1e-6, resolution=nothing,
                 kwargs...)
    if specaxis == :λ
        specx, Iω, zactual = getIω(output, specaxis, zslice, specrange=λrange, resolution=resolution)
    else
        specx, Iω, zactual = getIω(output, specaxis, zslice, resolution=resolution)
    end
    speclims, speclabel, specxfac = getspeclims(λrange, specaxis)
    multimode, modestrs = get_modes(output)
    if multimode
        modes = isnothing(modes) ? (1:size(Iω, 2)) : modes
        if modes == :sum
            Iω = dropdims(sum(Iω, dims=2), dims=2)
            modestrs = join(modestrs, "+")
            nmodes = 1
        else
            isnothing(modes) && (modes = 1:length(modestrs))
            Iω = Iω[:, modes, :]
            modestrs = modestrs[modes]
            nmodes = length(modes)
        end
    end

    specx .*= specxfac

    sfig = plt.figure()
    if multimode && nmodes > 1
        _plot_slice_mm(plt.gca(), specx, Iω, zactual, modestrs, log10; kwargs...)
    else
        zs = [@sprintf("%.2f cm", zi*100) for zi in zactual]
        label = multimode ? zs.*" ($modestrs)" : zs
        for iz in eachindex(zactual)
            (log10 ? plt.semilogy : plt.plot)(specx, Iω[:, iz]; label=label[iz], kwargs...)
        end
    end
    plt.legend(frameon=false)
    plt.xlabel(speclabel)
    plt.ylabel("Spectral energy density")
    log10 && plt.ylim(3*maximum(Iω)*log10min, 3*maximum(Iω))
    plt.xlim(speclims...)
    sfig.set_size_inches(8.5, 5)
    sfig.tight_layout()
    sfig
end

dashes = [(0, (10, 1)),
          (0, (5, 1)),
          (0, (1, 0.5)),
          (0, (1, 0.5, 1, 0.5, 3, 1)),
          (0, (5, 1, 1, 1))]

function _plot_slice_mm(ax, x, y, z, modestrs, log10=false, fwhm=false; kwargs...)
    pfun = (log10 ? ax.semilogy : ax.plot)
    for sidx = 1:size(y, 3) # iterate over z-slices
        zs = @sprintf("%.2f cm", z[sidx]*100)
        line = pfun(x, y[:, 1, sidx]; label="$zs ($(modestrs[1]))", kwargs...)[1]
        for midx = 2:size(y, 2) # iterate over modes
            pfun(x, y[:, midx, sidx], linestyle=dashes[midx], color=line.get_color(),
                 label="$zs ($(modestrs[midx]))"; kwargs...)
        end
    end
end

spectrogram(output::AbstractOutput, args...; kwargs...) = spectrogram(
    makegrid(output), output, args...; kwargs...)

function spectrogram(grid::Grid.AbstractGrid, Eω::AbstractArray, specaxis=:λ;
                     propagate=nothing, kwargs...)
    t, Et = getEt(grid, Eω; propagate=propagate, oversampling=1)
    spectrogram(t, Et, specaxis; kwargs...)
end

function spectrogram(grid::Grid.AbstractGrid, output, zslice, specaxis=:λ;
                     propagate=nothing, kwargs...)
    t, Et, zactual = getEt(output, zslice; oversampling=1, propagate=propagate)
    Et = Et[:, 1]
    spectrogram(t, Et, specaxis; kwargs...)
end

function spectrogram(t::AbstractArray, Et::AbstractArray, specaxis=:λ;
                     trange, N, fw, λrange=(150e-9, 2000e-9), log=false, dBmin=-40,
                     kwargs...)
    ω = Maths.rfftfreq(t)[2:end]
    tmin, tmax = extrema(trange)
    tg = collect(range(tmin, tmax, length=N))
    g = Maths.gabor(t, real(Et), tg, fw)
    g = g[2:end, :]

    specy, Ig = getIω(ω, g*Maths.rfftnorm(t[2]-t[1]), specaxis)
    speclims, speclabel, specyfac = getspeclims(λrange, specaxis)

    log && (Ig = 10*log10.(Maths.normbymax(Ig)))

    fig = plt.figure()
    plt.pcolormesh(tg.*1e15, specyfac*specy, Ig; shading="auto", kwargs...)
    plt.ylim(speclims...)
    plt.ylabel(speclabel)
    plt.xlabel("Time (fs)")
    log && plt.clim(dBmin, 0)
    plt.colorbar()
    fig
end

function energy(output; modes=nothing, bandpass=nothing, figsize=(7, 5))
    e = Processing.energy(output; bandpass=bandpass)
    eall = Processing.energy(output)

    multimode, modestrs = get_modes(output)
    if multimode
        e0 = sum(eall[:, 1])
        modes = isnothing(modes) ? (1:size(e, 1)) : modes
        if modes == :sum
            e = dropdims(sum(e, dims=1), dims=1)
            modestrs = join(modestrs, "+")
            nmodes = 1
        else
            isnothing(modes) && (modes = 1:length(modestrs))
            e = e[modes, :]
            modestrs = modestrs[modes]
            nmodes = length(modes)
        end
    else
        e0 = eall[1]
    end

    z = output["z"]*100

    fig = plt.figure()
    ax = plt.axes()
    ax.plot(z, 1e6*e')
    ax.set_xlim(extrema(z)...)
    ax.set_ylim(ymin=0)
    ax.set_xlabel("Distance (cm)")
    ax.set_ylabel("Energy (μJ)")
    rax = ax.twinx()
    rax.plot(z, 100*(e/e0)', linewidth=0)
    lims = ax.get_ylim()
    rax.set_ylim(100/(1e6*e0).*lims)
    rax.set_ylabel("Conversion efficiency (%)")
    fig.set_size_inches(figsize...)
    fig
end

function auto_fwhm_arrows(ax, x, y; color="k", arrowlength=nothing, hpad=0, linewidth=1,
                                    text=nothing, units="fs")
    left, right = Maths.level_xings(x, y)
    fw = abs(right - left)
    halfmax = maximum(y)/2
    arrowlength = isnothing(arrowlength) ? 2*fw : arrowlength

    ax.annotate("", xy=(left-hpad, halfmax),
                xytext=(left-hpad-arrowlength, halfmax),
                arrowprops=Dict("arrowstyle" => "->",
                                "color" => color,
                                "linewidth" => linewidth))
    ax.annotate("", xy=(right+hpad, halfmax),
                xytext=(right+hpad+arrowlength, halfmax),
                arrowprops=Dict("arrowstyle" => "->",
                                "color" => color,
                                "linewidth" => linewidth))

    if text == :left
        ax.text(left-arrowlength/2, 1.1*halfmax, @sprintf("%.1f %s", fw, units),
                ha="right", color=color)
    elseif text == :right
        ax.text(right+arrowlength/2, 1.1*halfmax, @sprintf("%.1f %s", fw, units),
                color=color)
    end
end

function add_fwhm_legends(ax, unit)
    leg = ax.get_legend()
    texts = leg.get_texts()
    handles, labels = ax.get_legend_handles_labels()

    for (ii, line) in enumerate(handles)
        xy = line.get_xydata()
        fw = Maths.fwhm(xy[:, 1], xy[:, 2])
        t = texts[ii]
        s = t.get_text()
        s *= @sprintf(" [%.2f %s]", fw, unit)
        t.set_text(s)
    end
end

end


text:
using Test;

using AdalmPluto;

@testset "libIIO/toplevel.jl" begin
    # disable the assertions
    toggleNoAssertions(true);

    # C_iio_library_get_version
     # @test C_iio_library_get_version() == (0, 21, "565bf68");

    # C_iio_get_backends_count
    @test C_iio_get_backends_count() > 0;

    # C_iio_has_backend
    @test C_iio_has_backend("usb") == true;
    @test C_iio_has_backend("dummy") == false;

    # C_iio_get_backend
    @test C_iio_get_backend(UInt32(0)) != "";
    @test C_iio_get_backend(UInt32(666)) == "";

    # C_iio_strerror
    @test C_iio_strerror(1) == "Operation not permitted (1)";
end

ext: jl


text:
using Turing: VarEstimator, add_sample!, get_var, reset!
using Test

let
    D = 1000
    ve = VarEstimator{Float64}(0, zeros(D), zeros(D))
    ve = reset!(add_sample!(ve, randn(D)))

    # Check that reseting zeros everything.
    @test ve.n === 0
    @test ve.μ == zeros(D)
    @test ve.M == zeros(D)

    # Ensure that asking for the variance doesn't mutate the VarEstimator.
    add_sample!(ve, randn(D))
    add_sample!(ve, randn(D))
    μ, M = deepcopy(ve.μ), deepcopy(ve.M)
    get_var(ve)
    @test ve.μ == μ
    @test ve.M == M
end

# Check that the estimated variance is approximately correct.
let
    D = 10
    ve = VarEstimator{Float64}(0, zeros(D), zeros(D))

    for _ = 1:10000
        s = randn(D)
        add_sample!(ve, s)
    end

    var = get_var(ve)

    @test var ≈ ones(D) atol=0.5
end


text:
struct PolyharmonicSpline
    dim::Int64
    order::Int64
    coeff::Vector{Float64}
    centers::Array{Float64,2}
    error::Float64
end

function polyharmonicK(r,K)
    if iseven(K)
        iszero(r) && return zero(r)
        if r >= 1.0
            return (r^K)*log(r)
        elseif 0.0 < r < 1.0
            return (r.^(K-1))*log(r.^r)
        elseif iszero(r) # Needed for autodiff to work at zero
            return zero(r)
        end
    else
        return r^K
    end
end

function PolyharmonicSpline(K::Int64, centers::Array{Float64,2}, values::Array{Float64}; s = 0.0)
    m,n = size(centers)
    m != length(values) && throw(DimensionMismatch())

    M = zeros(m,m)
    N = zeros(m,n+1)

    for i=1:m
        N[i,1] = 1
        N[i,2:end] = centers[i,:]
        for j=1:m
            M[i,j] = polyharmonicK(norm(centers[i,:] .- centers[j,:]),K)
        end
    end
    M = M + s*I
    L = vcat(hcat(M,N),hcat(N', zeros(n+1,n+1)))

    w = pinv(L)*vcat(values,zeros(n+1))

    ivalues = zeros(m)
    for i=1:m
        tmp = 0.0
        for j=1:m
            tmp = tmp + w[j]*polyharmonicK(norm(centers[i,:] .- centers[j,:]),K)
        end
        tmp = tmp + w[m+1]
        for j=2:n+1
            tmp = tmp + w[m+j]*centers[i,j-1]
        end
        ivalues[i] = tmp
    end
    error = norm(values .- ivalues)

    return PolyharmonicSpline(n,K,w,centers,error)
end

function PolyharmonicSpline(K::Int64, centers::Vector{Float64},values::Vector{Float64};s = 0.0)
    PolyharmonicSpline(K,reshape(centers,length(centers),1),values,s=s)
end

function (S::PolyharmonicSpline)(x::T...) where T <: Real
    n = length(x)
    n != S.dim && throw(DimensionMismatch("$n != $(S.dim)"))

    v = 0.0
    l = length(S.coeff)-(n+1)
    for j=1:l
        v = v + S.coeff[j]*polyharmonicK(norm(x .- S.centers[j,:]), S.order)
    end

    v = v + S.coeff[l+1]
    for j=2:n+1
        v = v + S.coeff[l+j]*x[j-1]
    end

    return v
end


text:
struct Wire <: StructuralElement
    a::NTuple{3,Float64}    # (x,y,z) of first end of wire (in meters)
    b::NTuple{3,Float64}    # (x,y,z) of second end of wire (in meters)
    rad::Float64            # radius of first wire segment (in meters)
    n_segments::Int         # number of segments

    # Advanced fields
    rdel::Float64   # ratio of length of segment to length of previous
    rrad::Float64   # ratio of radii of adjacent segments
end
Wire(a, b, rad, n_segments) = Wire(a, b, rad, n_segments, 1, 1)

# If there are 21 segments, they are indexed 1-21 (not 0-20).
#  So we'd want to return 11, and div(21,2) returns 10.
Base.middle(wire::Wire) = div(wire.n_segments, 2) + 1

function Base.length(w::Wire)
    dx = w.a[1] - w.b[1]
    dy = w.a[2] - w.b[2]
    dz = w.a[3] - w.b[3]

    return sqrt(dx*dx + dy*dy + dz*dz)
end


text:
module Day10

using Statistics

import Base.show, Base.-, Base.*
import ..data_dir # from parent module
input = read(joinpath(data_dir, "day10"), String)

export part1, part2

struct Bracket
    left::AbstractString
    right::AbstractString
    orientation::AbstractString
    function Bracket(left::AbstractString, right::AbstractString, orientation::AbstractString)
        if orientation ∉ ["left", "right"]
            throw(ArgumentError("orientation must be \"left\" or \"right\", got: $orientation"))
        end
        new(left, right, orientation)
    end
end

function Bracket(b::AbstractString)
    if b == "("
        Bracket("(", ")", "left")
    elseif b == ")"
        Bracket("(", ")", "right")
    elseif b == "["
        Bracket("[", "]", "left")
    elseif b == "]"
        Bracket("[", "]", "right")
    elseif b == "{"
        Bracket("{", "}", "left")
    elseif b == "}"
        Bracket("{", "}", "right")
    elseif b == "<"
        Bracket("<", ">", "left")
    elseif b == ">"
        Bracket("<", ">", "right")
    else
        throw(ArgumentError("No shortcut for bracket: $b"))
    end
end

function -(bracket::Bracket)
    if bracket.orientation == "left"
        return (Bracket(bracket.left, bracket.right, "right"))
    end
    Bracket(bracket.left, bracket.right, "left")
end

function strip_matching_brackets(brackets::Vector{Bracket})
    for i = 1:length(brackets)-1
        if brackets[i].orientation == "left" && brackets[i] == -brackets[i+1]
            stripped = [brackets[j] for j in 1:length(brackets) if j ∉ [i, i + 1]]
            return strip_matching_brackets(stripped)
        end
    end
    return brackets
end

# Input parsing
function convert_input_to_lines(input::AbstractString)
    raw_lines = split(input, "\n")
    [Bracket.(split(line, "")) for line in raw_lines]
end

# Part 1
const bad_bracket_scores = Dict(
    Bracket(")") => 3,
    Bracket("]") => 57,
    Bracket("}") => 1197,
    Bracket(">") => 25137
)

function first_bad_bracket(brackets::Vector{Bracket})
    stripped = strip_matching_brackets(brackets)
    for i = 2:length(stripped)
        if stripped[i].orientation == "right" && -stripped[i] != stripped[i-1]
            return (stripped[i])
        end
    end
    nothing
end

function part1(input = input)
    lines = convert_input_to_lines(input)
    first_bad_brackets = first_bad_bracket.(lines)
    sum([bad_bracket_scores[bracket] for bracket in first_bad_brackets if !isnothing(bracket)])
end

# Part 2
const completion_values = Dict(
    Bracket(")") => 1,
    Bracket("]") => 2,
    Bracket("}") => 3,
    Bracket(">") => 4
)

function calculate_completion_score(brackets::Vector{Bracket})
    stripped = strip_matching_brackets(brackets)
    if !isnothing(first_bad_bracket(stripped))
        return nothing
    end
    right_half = [-bracket for bracket in reverse(stripped)]
    completion_score = 0
    for bracket in right_half
        completion_score = completion_score * 5
        completion_score += completion_values[bracket]
    end
    completion_score
end

function part2(input = input)
    lines = convert_input_to_lines(input)
    completion_scores = calculate_completion_score.(lines)
    non_missing = [score for score in completion_scores if !isnothing(score)]
    Int(median(non_missing))
end

end #module


text:
module VariationalAE

export VAE, fit!, predict, generate, reconstruct

include("core.jl")

end # module


text:
## FUNCTIONS TO LOAD MODEL IMPLEMENTATION CODE

### Helpers

function _append!(program, ex, doprint::Bool, tick_early::Bool)
    str = string(ex)
    doprint && push!(program.args, :(print($str)))
    doprint && !tick_early && push!(program.args, :(println(" \u2714")))
    push!(program.args, ex)
    tick_early && doprint && push!(program.args, :(println(" \u2714")))
    return program
end

"""
    load_path(model::String, pkg=nothing)

Return the load path for model type with name `model`, specifying the
package name `pkg` to resolve name conflicts if necessary.

    load_path(model)

Return the load path of a `model` instance or type. Usually requires
necessary model code to have been separately loaded. Supply a string
as above if code is not loaded.

"""
function MLJModelInterface.load_path(proxy::ModelProxy)
    handle = (name=proxy.name, pkg=proxy.package_name)
    return INFO_GIVEN_HANDLE[handle][:load_path]
end
function MLJModelInterface.load_path(name::String; pkg=nothing)
    proxy = info(name; pkg=pkg)
    return load_path(proxy)
end

# to also get pkg, which could be different from glue code pkg
# appearing in load_path:
function load_path_and_pkg(name::String; pkg=nothing)
    proxy = info(name; pkg=pkg)
    handle = (name=proxy.name, pkg=proxy.package_name)
    _info = INFO_GIVEN_HANDLE[handle]

    return _info[:load_path], _info[:package_name]
end

"""
    @load name pkg=nothing verbosity=0 name=nothing scope=:global install_pkgs=false

Load the model implementation code for the model named in the first
argument into the calling module, specfying `pkg` in the case of
duplicate names. Return a model instance with default
hyper-parameters.

To bind the new model type to a name different from the first
argument, specify `name=...`.

### Examples

    @load DecisionTreeeRegressor
    @load PCA verbosity=1
    @load SVC pkg=LIBSVM name=MyOtherSVC

See also [`load`](@ref)

"""
macro load(name_ex, kw_exs...)
    program, instance_prgm = _load(__module__, name_ex, kw_exs...)
    append!(program.args, instance_prgm.args)
    esc(program)
end

macro loadcode(name_ex, kw_exs...)
    program, _ = _load(__module__, name_ex, kw_exs...)
    esc(program)
end

# builds the program to be evaluated by the @load macro:
function _load(modl, name_ex, kw_exs...)

    # initialize:
    program = quote end
    instance_prgm = quote end

    # fallbacks:
    pkg = nothing
    verbosity = 1
    new_name = nothing
    scope = :global
    install_pkgs = false

    # parse name_ex:
    name_ = string(name_ex)
    # get rid of parentheses in `name_`, as in
    # "(MLJModels.Clustering).KMedoids":
    name = filter(name_) do c !(c in ['(',')']) end

    # parse kwargs:
    warning = "Invalid @load syntax.\n "*
    "Sample usage: @load PCA pkg=\"MultivariateStats\" verbosity=1"
    for ex in kw_exs
        ex.head == :(=) || throw(ArgumentError(warning))
        variable_ex = ex.args[1]
        value_ex = ex.args[2]
        if variable_ex == :pkg
            pkg = string(value_ex)
        elseif variable_ex == :verbosity
            verbosity = value_ex
        elseif variable_ex == :name
            new_name = value_ex
        elseif variable_ex == :scope
            scope = value_ex
        elseif variable_ex == :install_pkgs
            install_pkgs = value_ex
        else
            throw(ArgumentError(warning))
        end
    end

    if scope == :(:global)
        scope = :global
    end
    if scope == :(:local)
        scope = :local
    end
    if !( scope in [:global, :local] )
        throw(ArgumentError("Invalid value for scope: $(scope). Valid values are :global or :local"))
    end

    # are we printing stuff to stdout?
    doprint = verbosity > 0

    doprint && @info "For silent loading, specify `verbosity=0`. "

    # get load path and update pkg (could be `nothing`):
    path, pkg = load_path_and_pkg(name, pkg=pkg)

    # see if the model type has already been loaded:
    handle = (name=name, pkg=pkg)
    type_already_loaded = handle in map(localmodels(modl=modl)) do m
        (name=m.name, pkg=m.package_name)
    end

    # if so, return with program generating an instance:
    if type_already_loaded && new_name == nothing
        doprint && @info "Model code for $name already loaded"
        # return an instance of the type:
        for M in localmodeltypes(modl)
            i = info(M)
            if i.name == name && i.package_name == pkg
                _append!(instance_prgm, :($M()), doprint, false)
                return program, instance_prgm
            end
        end
    end

    # determine `new_name`, to be bound to imported model type (in
    # general, different from `name`):
    if new_name === nothing
        new_name = MLJBase.available_name(modl, Symbol(name))
        new_name == Symbol(name) || verbosity < 0 ||
            @warn "New model type being bound to "*
            "`$new_name` to avoid conflict with an existing name. "
    else
        new_name = Symbol(new_name)
    end

    path_components = Symbol.(split(path, '.') )

    # get pkg containing implementation of model API implementation:
    api_pkg = path_components[1]
    pkg = Symbol(pkg)

    # if needed, put MLJModels in the calling module's namespace:
    if api_pkg == :MLJModels
        load_ex =
            isdefined(modl, :MLJ) ? :(import MLJ.MLJModels) :
            :(import MLJModels)
        _append!(program, load_ex, doprint, true)
        # TODO: remove next line of code after disintegration of
        # MLJModels (for triggering loading of glue code module):
        api_pkg == pkg || _append!(program, :(import $pkg), doprint, true)
    end

    import_ex = :(import $api_pkg)
    if install_pkgs
        install_ex = quote
            try
                $(import_ex)
            catch
                import Pkg
                Pkg.add($(string(api_pkg)))
            end
        end
        install_ex = MacroTools.prewalk(rmlines, install_ex)
        _append!(program, install_ex, doprint, true)
    end

    path_ex = path |> Meta.parse
    api_pkg == :MLJmodels || _append!(program, import_ex, doprint, true)
    if scope == :local
        # we cannot use the `const` keyword in front of local variables
        _append!(program, :($new_name = $path_ex), doprint, true)
    else
        _append!(program, :(const $new_name = $path_ex), doprint, true)
    end

    instance_ex = doprint ? :($new_name()) : :($new_name();)
    _append!(instance_prgm, instance_ex, doprint, false)

    return program, instance_prgm
end

## NO LONGER SUPPORTED

_deperror() = error(
    "The `load` function is no longer supported. "*
    "Use the `@load` macro instead, as in "*
    "`@load RandomForestRegressor pkg = DecisionTree`.\n"*
    "For explicit importing, you can discover a model's "*
    "full load path with the `load_path` function, as in "*
    "`load_path(\"RandomForestRegressor\", pkg=\"DecisionTree\")`. )")

load(proxy::ModelProxy; kwargs...) = _deperror()
load(name::String; kwargs...) = _deperror()


text:
@doc raw"""
    BinaryRBM(a, b, w)
    BinaryRBM(N, M)

Construct an RBM with binary visible and hidden units, which has an energy function:

```math
E(v, h) = -a'v - b'h - v'wh
```

Equivalent to `RBM(Binary(a), Binary(b), w)`.
"""
function BinaryRBM(a::AbstractArray, b::AbstractArray, w::AbstractArray)
    @assert size(w) == (size(a)..., size(b)...)
    return RBM(Binary(a), Binary(b), w)
end

function BinaryRBM(::Type{T}, N::Union{Int,TupleN{Int}}, M::Union{Int,TupleN{Int}}) where {T}
    a = zeros(T, N...)
    b = zeros(T, M...)
    w = zeros(T, N..., M...)
    return BinaryRBM(a, b, w)
end

BinaryRBM(N::Union{Int,TupleN{Int}}, M::Union{Int,Tuple{Vararg{Int}}}) = BinaryRBM(Float64, N, M)

Base.repeat(l::Binary, n::Int...) = Binary(repeat(l.θ, n...))


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: ivs
using AWS.Compat
using AWS.UUIDs

"""
    batch_get_channel(arns)
    batch_get_channel(arns, params::Dict{String,<:Any})

Performs GetChannel on multiple ARNs simultaneously.

# Arguments
- `arns`: Array of ARNs, one per channel.

"""
batch_get_channel(arns; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/BatchGetChannel", Dict{String, Any}("arns"=>arns); aws_config=aws_config)
batch_get_channel(arns, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/BatchGetChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arns"=>arns), params)); aws_config=aws_config)

"""
    batch_get_stream_key(arns)
    batch_get_stream_key(arns, params::Dict{String,<:Any})

Performs GetStreamKey on multiple ARNs simultaneously.

# Arguments
- `arns`: Array of ARNs, one per channel.

"""
batch_get_stream_key(arns; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/BatchGetStreamKey", Dict{String, Any}("arns"=>arns); aws_config=aws_config)
batch_get_stream_key(arns, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/BatchGetStreamKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arns"=>arns), params)); aws_config=aws_config)

"""
    create_channel()
    create_channel(params::Dict{String,<:Any})

Creates a new channel and an associated stream key to start streaming.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authorized"`: Whether the channel is private (enabled for playback authorization).
  Default: false.
- `"latencyMode"`: Channel latency mode. Use NORMAL to broadcast and deliver live video up
  to Full HD. Use LOW for near-real-time interaction with viewers. (Note: In the Amazon IVS
  console, LOW and NORMAL correspond to Ultra-low and Standard, respectively.) Default: LOW.
- `"name"`: Channel name.
- `"recordingConfigurationArn"`: Recording-configuration ARN. Default: \"\" (empty string,
  recording is disabled).
- `"tags"`: Array of 1-50 maps, each of the form string:string (key:value).
- `"type"`: Channel type, which determines the allowable resolution and bitrate. If you
  exceed the allowable resolution or bitrate, the stream probably will disconnect
  immediately. Default: STANDARD. Valid values:    STANDARD: Multiple qualities are generated
  from the original input, to automatically give viewers the best experience for their
  devices and network conditions. Vertical resolution can be up to 1080 and bitrate can be up
  to 8.5 Mbps.    BASIC: Amazon IVS delivers the original input to viewers. The viewer’s
  video-quality choice is limited to the original input. Vertical resolution can be up to 480
  and bitrate can be up to 1.5 Mbps.
"""
create_channel(; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/CreateChannel"; aws_config=aws_config)
create_channel(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/CreateChannel", params; aws_config=aws_config)

"""
    create_recording_configuration(destination_configuration)
    create_recording_configuration(destination_configuration, params::Dict{String,<:Any})

Creates a new recording configuration, used to enable recording to Amazon S3.  Known issue:
In the us-east-1 region, if you use the AWS CLI to create a recording configuration, it
returns success even if the S3 bucket is in a different region. In this case, the state of
the recording configuration is CREATE_FAILED (instead of ACTIVE). (In other regions, the
CLI correctly returns failure if the bucket is in a different region.)  Workaround: Ensure
that your S3 bucket is in the same region as the recording configuration. If you create a
recording configuration in a different region as your S3 bucket, delete that recording
configuration and create a new one with an S3 bucket from the correct region.

# Arguments
- `destination_configuration`: A complex type that contains a destination configuration for
  where recorded video will be stored.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"name"`: An arbitrary string (a nickname) that helps the customer identify that
  resource. The value does not need to be unique.
- `"tags"`: Array of 1-50 maps, each of the form string:string (key:value).
"""
create_recording_configuration(destinationConfiguration; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/CreateRecordingConfiguration", Dict{String, Any}("destinationConfiguration"=>destinationConfiguration); aws_config=aws_config)
create_recording_configuration(destinationConfiguration, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/CreateRecordingConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("destinationConfiguration"=>destinationConfiguration), params)); aws_config=aws_config)

"""
    create_stream_key(channel_arn)
    create_stream_key(channel_arn, params::Dict{String,<:Any})

Creates a stream key, used to initiate a stream, for the specified channel ARN. Note that
CreateChannel creates a stream key. If you subsequently use CreateStreamKey on the same
channel, it will fail because a stream key already exists and there is a limit of 1 stream
key per channel. To reset the stream key on a channel, use DeleteStreamKey and then
CreateStreamKey.

# Arguments
- `channel_arn`: ARN of the channel for which to create the stream key.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tags"`: Array of 1-50 maps, each of the form string:string (key:value).
"""
create_stream_key(channelArn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/CreateStreamKey", Dict{String, Any}("channelArn"=>channelArn); aws_config=aws_config)
create_stream_key(channelArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/CreateStreamKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn), params)); aws_config=aws_config)

"""
    delete_channel(arn)
    delete_channel(arn, params::Dict{String,<:Any})

Deletes the specified channel and its associated stream keys. If you try to delete a live
channel, you will get an error (409 ConflictException). To delete a channel that is live,
call StopStream, wait for the Amazon EventBridge \"Stream End\" event (to verify that the
stream's state was changed from Live to Offline), then call DeleteChannel. (See  Using
EventBridge with Amazon IVS.)

# Arguments
- `arn`: ARN of the channel to be deleted.

"""
delete_channel(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/DeleteChannel", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_channel(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/DeleteChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    delete_playback_key_pair(arn)
    delete_playback_key_pair(arn, params::Dict{String,<:Any})

Deletes a specified authorization key pair. This invalidates future viewer tokens generated
using the key pair’s privateKey. For more information, see Setting Up Private Channels in
the Amazon IVS User Guide.

# Arguments
- `arn`: ARN of the key pair to be deleted.

"""
delete_playback_key_pair(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/DeletePlaybackKeyPair", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_playback_key_pair(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/DeletePlaybackKeyPair", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    delete_recording_configuration(arn)
    delete_recording_configuration(arn, params::Dict{String,<:Any})

Deletes the recording configuration for the specified ARN. If you try to delete a recording
configuration that is associated with a channel, you will get an error (409
ConflictException). To avoid this, for all channels that reference the recording
configuration, first use UpdateChannel to set the recordingConfigurationArn field to an
empty string, then use DeleteRecordingConfiguration.

# Arguments
- `arn`: ARN of the recording configuration to be deleted.

"""
delete_recording_configuration(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/DeleteRecordingConfiguration", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_recording_configuration(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/DeleteRecordingConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    delete_stream_key(arn)
    delete_stream_key(arn, params::Dict{String,<:Any})

Deletes the stream key for the specified ARN, so it can no longer be used to stream.

# Arguments
- `arn`: ARN of the stream key to be deleted.

"""
delete_stream_key(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/DeleteStreamKey", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
delete_stream_key(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/DeleteStreamKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    get_channel(arn)
    get_channel(arn, params::Dict{String,<:Any})

Gets the channel configuration for the specified channel ARN. See also BatchGetChannel.

# Arguments
- `arn`: ARN of the channel for which the configuration is to be retrieved.

"""
get_channel(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetChannel", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_channel(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    get_playback_key_pair(arn)
    get_playback_key_pair(arn, params::Dict{String,<:Any})

Gets a specified playback authorization key pair and returns the arn and fingerprint. The
privateKey held by the caller can be used to generate viewer authorization tokens, to grant
viewers access to private channels. For more information, see Setting Up Private Channels
in the Amazon IVS User Guide.

# Arguments
- `arn`: ARN of the key pair to be returned.

"""
get_playback_key_pair(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetPlaybackKeyPair", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_playback_key_pair(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetPlaybackKeyPair", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    get_recording_configuration(arn)
    get_recording_configuration(arn, params::Dict{String,<:Any})

Gets the recording configuration for the specified ARN.

# Arguments
- `arn`: ARN of the recording configuration to be retrieved.

"""
get_recording_configuration(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetRecordingConfiguration", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_recording_configuration(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetRecordingConfiguration", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    get_stream(channel_arn)
    get_stream(channel_arn, params::Dict{String,<:Any})

Gets information about the active (live) stream on a specified channel.

# Arguments
- `channel_arn`: Channel ARN for stream to be accessed.

"""
get_stream(channelArn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetStream", Dict{String, Any}("channelArn"=>channelArn); aws_config=aws_config)
get_stream(channelArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetStream", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn), params)); aws_config=aws_config)

"""
    get_stream_key(arn)
    get_stream_key(arn, params::Dict{String,<:Any})

Gets stream-key information for a specified ARN.

# Arguments
- `arn`: ARN for the stream key to be retrieved.

"""
get_stream_key(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetStreamKey", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
get_stream_key(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/GetStreamKey", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)

"""
    import_playback_key_pair(public_key_material)
    import_playback_key_pair(public_key_material, params::Dict{String,<:Any})

Imports the public portion of a new key pair and returns its arn and fingerprint. The
privateKey can then be used to generate viewer authorization tokens, to grant viewers
access to private channels. For more information, see Setting Up Private Channels in the
Amazon IVS User Guide.

# Arguments
- `public_key_material`: The public portion of a customer-generated key pair.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"name"`: An arbitrary string (a nickname) assigned to a playback key pair that helps the
  customer identify that resource. The value does not need to be unique.
- `"tags"`: Any tags provided with the request are added to the playback key pair tags.
"""
import_playback_key_pair(publicKeyMaterial; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ImportPlaybackKeyPair", Dict{String, Any}("publicKeyMaterial"=>publicKeyMaterial); aws_config=aws_config)
import_playback_key_pair(publicKeyMaterial, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ImportPlaybackKeyPair", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("publicKeyMaterial"=>publicKeyMaterial), params)); aws_config=aws_config)

"""
    list_channels()
    list_channels(params::Dict{String,<:Any})

Gets summary information about all channels in your account, in the AWS region where the
API request is processed. This list can be filtered to match a specified name or
recording-configuration ARN. Filters are mutually exclusive and cannot be used together. If
you try to use both filters, you will get an error (409 ConflictException).

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"filterByName"`: Filters the channel list to match the specified name.
- `"filterByRecordingConfigurationArn"`: Filters the channel list to match the specified
  recording-configuration ARN.
- `"maxResults"`: Maximum number of channels to return. Default: 50.
- `"nextToken"`: The first channel to retrieve. This is used for pagination; see the
  nextToken response field.
"""
list_channels(; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListChannels"; aws_config=aws_config)
list_channels(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListChannels", params; aws_config=aws_config)

"""
    list_playback_key_pairs()
    list_playback_key_pairs(params::Dict{String,<:Any})

Gets summary information about playback key pairs. For more information, see Setting Up
Private Channels in the Amazon IVS User Guide.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The first key pair to retrieve. This is used for pagination; see the
  nextToken response field. Default: 50.
- `"nextToken"`: Maximum number of key pairs to return.
"""
list_playback_key_pairs(; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListPlaybackKeyPairs"; aws_config=aws_config)
list_playback_key_pairs(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListPlaybackKeyPairs", params; aws_config=aws_config)

"""
    list_recording_configurations()
    list_recording_configurations(params::Dict{String,<:Any})

Gets summary information about all recording configurations in your account, in the AWS
region where the API request is processed.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of recording configurations to return. Default: 50.
- `"nextToken"`: The first recording configuration to retrieve. This is used for
  pagination; see the nextToken response field.
"""
list_recording_configurations(; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListRecordingConfigurations"; aws_config=aws_config)
list_recording_configurations(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListRecordingConfigurations", params; aws_config=aws_config)

"""
    list_stream_keys(channel_arn)
    list_stream_keys(channel_arn, params::Dict{String,<:Any})

Gets summary information about stream keys for the specified channel.

# Arguments
- `channel_arn`: Channel ARN used to filter the list.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of streamKeys to return. Default: 50.
- `"nextToken"`: The first stream key to retrieve. This is used for pagination; see the
  nextToken response field.
"""
list_stream_keys(channelArn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListStreamKeys", Dict{String, Any}("channelArn"=>channelArn); aws_config=aws_config)
list_stream_keys(channelArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListStreamKeys", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn), params)); aws_config=aws_config)

"""
    list_streams()
    list_streams(params::Dict{String,<:Any})

Gets summary information about live streams in your account, in the AWS region where the
API request is processed.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of streams to return. Default: 50.
- `"nextToken"`: The first stream to retrieve. This is used for pagination; see the
  nextToken response field.
"""
list_streams(; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListStreams"; aws_config=aws_config)
list_streams(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/ListStreams", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Gets information about AWS tags for the specified ARN.

# Arguments
- `resource_arn`: The ARN of the resource to be retrieved.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of tags to return. Default: 50.
- `"nextToken"`: The first tag to retrieve. This is used for pagination; see the nextToken
  response field.
"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    put_metadata(channel_arn, metadata)
    put_metadata(channel_arn, metadata, params::Dict{String,<:Any})

Inserts metadata into the active stream of the specified channel. A maximum of 5 requests
per second per channel is allowed, each with a maximum 1 KB payload. (If 5 TPS is not
sufficient for your needs, we recommend batching your data into a single PutMetadata call.)
Also see Embedding Metadata within a Video Stream in the Amazon IVS User Guide.

# Arguments
- `channel_arn`: ARN of the channel into which metadata is inserted. This channel must have
  an active stream.
- `metadata`: Metadata to insert into the stream. Maximum: 1 KB per request.

"""
put_metadata(channelArn, metadata; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/PutMetadata", Dict{String, Any}("channelArn"=>channelArn, "metadata"=>metadata); aws_config=aws_config)
put_metadata(channelArn, metadata, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/PutMetadata", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn, "metadata"=>metadata), params)); aws_config=aws_config)

"""
    stop_stream(channel_arn)
    stop_stream(channel_arn, params::Dict{String,<:Any})

Disconnects the incoming RTMPS stream for the specified channel. Can be used in conjunction
with DeleteStreamKey to prevent further streaming to a channel.  Many streaming
client-software libraries automatically reconnect a dropped RTMPS session, so to stop the
stream permanently, you may want to first revoke the streamKey attached to the channel.

# Arguments
- `channel_arn`: ARN of the channel for which the stream is to be stopped.

"""
stop_stream(channelArn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/StopStream", Dict{String, Any}("channelArn"=>channelArn); aws_config=aws_config)
stop_stream(channelArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/StopStream", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("channelArn"=>channelArn), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds or updates tags for the AWS resource with the specified ARN.

# Arguments
- `resource_arn`: ARN of the resource for which tags are to be added or updated.
- `tags`: Array of tags to be added or updated.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes tags from the resource with the specified ARN.

# Arguments
- `resource_arn`: ARN of the resource for which tags are to be removed.
- `tag_keys`: Array of tags to be removed.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_channel(arn)
    update_channel(arn, params::Dict{String,<:Any})

Updates a channel's configuration. This does not affect an ongoing stream of this channel.
You must stop and restart the stream for the changes to take effect.

# Arguments
- `arn`: ARN of the channel to be updated.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authorized"`: Whether the channel is private (enabled for playback authorization).
- `"latencyMode"`: Channel latency mode. Use NORMAL to broadcast and deliver live video up
  to Full HD. Use LOW for near-real-time interaction with viewers. (Note: In the Amazon IVS
  console, LOW and NORMAL correspond to Ultra-low and Standard, respectively.)
- `"name"`: Channel name.
- `"recordingConfigurationArn"`: Recording-configuration ARN. If this is set to an empty
  string, recording is disabled. A value other than an empty string indicates that recording
  is enabled
- `"type"`: Channel type, which determines the allowable resolution and bitrate. If you
  exceed the allowable resolution or bitrate, the stream probably will disconnect
  immediately. Valid values:    STANDARD: Multiple qualities are generated from the original
  input, to automatically give viewers the best experience for their devices and network
  conditions. Vertical resolution can be up to 1080 and bitrate can be up to 8.5 Mbps.
  BASIC: Amazon IVS delivers the original input to viewers. The viewer’s video-quality
  choice is limited to the original input. Vertical resolution can be up to 480 and bitrate
  can be up to 1.5 Mbps.
"""
update_channel(arn; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/UpdateChannel", Dict{String, Any}("arn"=>arn); aws_config=aws_config)
update_channel(arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = ivs("POST", "/UpdateChannel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("arn"=>arn), params)); aws_config=aws_config)


text:
#=
Given an array of strictly the characters 'R', 'G', and 'B', segregate the values of the array so that all the Rs come first, the Gs come second, and the Bs come last. You can only swap elements of the array.

Do this in linear time and in-place.

For example, given the array ['G', 'B', 'R', 'R', 'B', 'R', 'G'], it should become ['R', 'R', 'R', 'G', 'G', 'B', 'B'].
=#

using Test
include("Solutions/problem35_segregateRGB.jl")

@test segregate_rgb(['G', 'B', 'R', 'R', 'B', 'R', 'G']) == ['R', 'R', 'R', 'G', 'G', 'B', 'B']

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# This file is a part of LegendDataTypes.jl, licensed under the MIT License (MIT).

export Geant4HDF5Input
struct Geant4HDF5Input <: AbstractLegendInput
    hdf5file::HDF5.HDF5File
end

struct GEARS_HDF5Input <: AbstractLegendInput
hdf5file::HDF5.HDF5File
end

struct G4SIMPLE_HDF5Input <: AbstractLegendInput
hdf5file::HDF5.HDF5File
end

Base.open(filename::AbstractString, ::Type{Geant4HDF5Input}) =
    Geant4HDF5Input(HDF5.h5open(filename, "r"))

Base.close(input::Geant4HDF5Input) = close(input.hdf5file)

function Base.getindex(input::Geant4HDF5Input, ::Colon)
    read(input)
end

function _d_open_if_exists(g, d_name::String)
    HDF5.exists(g, d_name) ? HDF5.d_open(g, d_name) : error("data set \"$d_name\" does not exist")
end

function Base.read(input::Geant4HDF5Input)
    if HDF5.exists(input.hdf5file, "/default_ntuples/t/")
        read(GEARS_HDF5Input(input.hdf5file))
    elseif HDF5.exists(input.hdf5file, "/default_ntuples/g4sntuple/")
        read(G4SIMPLE_HDF5Input(input.hdf5file))
    end
end

function Base.read(input::GEARS_HDF5Input)

    function _get_indices_raw(g)
        len = filter(x->x!=0x00, _d_open_if_exists(g, "e/pages")[:])
        indices_raw = []
        start = 0
        for i in len
            start+=1
            stop = start+i-1
            push!(indices_raw, collect(start:stop))
            start = stop
        end
        indices_raw
    end

    h5f = input.hdf5file
    g = HDF5.g_open(h5f, "default_ntuples/t/")
    indices_raw = _get_indices_raw(g)
    h5_edep = _d_open_if_exists(g, "e_data/pages")[:]
    indices_non_zero_energies  = filter(i-> sum(h5_edep[i]) > 0.0, indices_raw) # Filters events that don't deposit any energy
    indices = vcat(indices_non_zero_energies...)
    n_ind = length(indices)

    evtno = Int32.(vcat([ [j for i in 1:length(indices_non_zero_energies[j])] for j in 1:length(indices_non_zero_energies) ]...))
    detno = try Int32.(_d_open_if_exists(g, "vlm_data/pages")[:][indices]) catch ; ones(Int32, n_ind) end
    thit = ( try Float32.(_d_open_if_exists(g, "t0_data/pages")[:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"s"
    edep = ( try Float32.(_d_open_if_exists(g, "e_data/pages")[:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"keV"

    x0 = try Float32.(_d_open_if_exists(g, "x0_data/pages")[:][indices]) catch ; zeros(Float32, n_ind) end
    y0 = try Float32.(_d_open_if_exists(g, "y0_data/pages")[:][indices]) catch ; zeros(Float32, n_ind) end
    z0 = try Float32.(_d_open_if_exists(g, "z0_data/pages")[:][indices]) catch ; zeros(Float32, n_ind) end

    pos = [ SVector{3}(([ x0[i], y0[i], z0[i] ] .* u"mm")...) for i in 1:n_ind ]

    ekin = ( try Float32.(_d_open_if_exists(g, "k_data/pages")[:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"keV"# kinetic Energy of track [keV]
    stp = try Int32.(_d_open_if_exists(g, "stp_data/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # step number
    l = ( try Float32.(_d_open_if_exists(g, "l_data/pages")[:][indices]) catch ; zeros(Float32, n_ind) end  ) .* u"mm"# step length
    mom = try Int32.(_d_open_if_exists(g, "mom_data/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # parent id
    trk = try Int32.(_d_open_if_exists(g, "trk_data/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # trk id
    pdg = try Int32.(_d_open_if_exists(g, "pdg_data/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # particle id
    pro = try Int32.(_d_open_if_exists(g, "pro_data/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # process id

    hits = TypedTables.Table(
        evtno = evtno,
        detno = detno,
        thit = thit,
        edep = edep,
        pos = pos,
        ekin = ekin,
        stp = stp,
        l = l,
        mom = mom,
        trk = trk,
        pdg = pdg,
        pro = pro,
    )
    return hits
end

function Base.read(input::G4SIMPLE_HDF5Input)
    function _get_indices_raw(g)
        len = filter(x->x!=0x00, _d_open_if_exists(g, "e/pages")[:])
        indices_raw = []
        start = 0
        for i in len
            start+=1
            stop = start+i-1
            push!(indices_raw, collect(start:stop))
            start = stop
        end
        indices_raw
    end

    h5f = input.hdf5file
    g = HDF5.g_open(h5f, "default_ntuples/g4sntuple/")
    # indices_raw = _get_indices_raw(g)
    h5_edep = _d_open_if_exists(g, "event/pages")[:]
    # indices_non_zero_energies  = filter(i-> sum(h5_edep[i]) > 0.0, indices_raw) # Filters events that don't deposit any energy
    # indices = vcat(indices_non_zero_energies...)
    n_ind = length(h5_edep)
    indices = [i for i in 1:n_ind]

    evtno = try Int32.(_d_open_if_exists(g, "event/pages")[:][indices]) catch ; ones(Int32, n_ind) end
    irep = try Int32.(_d_open_if_exists(g, "iRep/pages")[:][indices]) catch ; ones(Int32, n_ind) end
    detno = try Int32.(_d_open_if_exists(g, "volID/pages")[:][indices]) catch ; ones(Int32, n_ind) end
    thit = ( try Float32.(_d_open_if_exists(g, "t/pages")[:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"s"
    edep = ( try Float32.(_d_open_if_exists(g, "Edep/pages")[:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"keV"

    x0 = try Float32.(_d_open_if_exists(g, "x/pages")[:][indices]) catch ; zeros(Float32, n_ind) end
    y0 = try Float32.(_d_open_if_exists(g, "y/pages")[:][indices]) catch ; zeros(Float32, n_ind) end
    z0 = try Float32.(_d_open_if_exists(g, "z/pages")[:][indices]) catch ; zeros(Float32, n_ind) end

    pos = [ SVector{3}(([ x0[i], y0[i], z0[i] ] .* u"mm")...) for i in 1:n_ind ]

    ekin = ( try Float32.(_d_open_if_exists(g, "KE_data/pages")[:][indices]) catch ; zeros(Float32, n_ind) end ) .* u"keV"# kinetic Energy of track [keV]
    stp = try Int32.(_d_open_if_exists(g, "step/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # step number
    mom = try Int32.(_d_open_if_exists(g, "parentID/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # parent id
    trk = try Int32.(_d_open_if_exists(g, "trackID/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # trk id
    pdg = try Int32.(_d_open_if_exists(g, "pid/pages")[:][indices]) catch ; zeros(Int32, n_ind) end # particle id

    hits = TypedTables.Table(
        evtno = evtno,
        detno = irep,
        thit = thit,
        edep = edep,
        pos = pos,
        ekin = ekin,
        volID = detno,
        stp = stp,
        mom = mom,
        trk = trk,
        pdg = pdg,
    )
    return hits
end


text:
struct TransmissionNetworkDistribution <: AbstractTN
  external::Array{Float64, 1}
  internal::Array{Float64, 2}

  function TransmissionNetworkDistribution(e::Vector{Float64},
                                           i::Array{Float64, 2})
    @boundscheck if !(length(e) == size(i, 1) == size(i, 2))
      error("Argument dimensions mismatched")
    end
    return new(e, i)
  end
end
function TransmissionNetworkDistribution(x::Vector{TransmissionNetwork})
  return @inbounds TransmissionNetworkDistribution(mean([y.external for y in x]), mean([y.internal for y in x]))
end

function TransmissionNetworkDistribution(x::Vector{TransmissionNetworkDistribution})
  return @inbounds TransmissionNetworkDistribution(mean([y.external for y in x]), mean([y.internal for y in x]))
end

const TNDistribution = TransmissionNetworkDistribution
const TransmissionNetworkPrior = TransmissionNetworkDistribution
const TNPrior = TransmissionNetworkPrior
const TransmissionNetworkPosterior = TransmissionNetworkDistribution
const TNPosterior = TransmissionNetworkPosterior

"""
Generate the posterior mode transmission network from a `TNDistribution`
"""
function mode(tnd::TNDistribution)
  tn = TransmissionNetwork(individuals(tnd))
  for i = 1:individuals(tnd)
    txfreq = [tnd.external[i]; tnd.internal[:,i]]
    if any(txfreq .> 0.)
      source = findmax([tnd.external[i]; tnd.internal[:,i]])[2]
      if source == 1
        tn.external[i] = 1
      else
        tn.internal[source-1, i] = 1
      end
    end
  end
  return tn
end


text:
using PersistenceDiagrams
using Ripserer
using RecipesBase
using Test

using Ripserer: plottable, ChainElement
using RecipesBase: apply_recipe

# Hack to avoid having to import Plots.
RecipesBase.is_key_supported(::Symbol) = true

"""
Idea: apply recipe and check the number of series on plots.
Not a perfect way to test, but at least it makes sure all points are plotted and checks that
there are no errors.
"""
series(args...; kwargs...) = apply_recipe(Dict{Symbol,Any}(kwargs), args...)

@testset "simplex plots" begin
    @testset "plottable" begin
        @test isequal(
            plottable(Simplex{0}(3, 1), 1:10), ([3], [:seriestype => :scatter], 0)
        )
        @test isequal(
            plottable(Simplex{1}(1, 1), 1:10), ([2, 1, NaN], [:seriestype => :path], 1)
        )
        @test isequal(
            plottable(Simplex{2}(3, 1), 1:10),
            ([4, 3, 1, 4, NaN], [:seriestype => :path], 2),
        )
        @test isequal(
            plottable(Simplex{3}(1, 1), 1:10),
            (
                [4, 3, 2, 4, 4, 3, 1, 4, 4, 2, 1, 4, 3, 2, 1, 3, NaN],
                [:seriestype => :path],
                3,
            ),
        )
    end
    @testset "plot" begin
        # No idea how to test this better. At least it checks there are no erros.
        data = collect(1:100)
        for dim in 0:3
            sx = Simplex{dim}(1, 1)
            fsx = Simplex{dim + 1}(1, 1)
            @test length(series(sx, data)) == 1
            @test length(series([sx], data)) == 1
            @test length(series([sx], data)) == 1
            @test length(
                series(
                    PersistenceInterval(
                        1.0,
                        1.0;
                        birth_simplex=sx,
                        death_simplex=fsx,
                        representative=[ChainElement{typeof(sx),typeof(1//1)}(sx, 1//1)],
                    ),
                    data,
                ),
            ) == 1
            @test_throws ErrorException series(PersistenceInterval(1.0, 1.0), data)
        end
    end
end


text:
using ApproxFun, Base.Test

l=u=1
λ=μ=1
N=M=10
cols=rows=1:N
data=ones(λ+μ+1,(l+u+1)*sum(cols))
A=ApproxFun.BandedBlockBandedMatrix(data,l,u,λ,μ,rows,cols)
v=ones(size(A,1))
M=full(A)
@test norm(A*v-M*v) ≤ 100eps()

A=ApproxFun.bbbrand(Float64,1,1,1,1,1:10,1:10)
B=ApproxFun.bbbrand(Float64,1,1,1,1,1:10,1:10)
@test_approx_eq A*B full(A)*full(B)

cols=Int[rand(1:k+2) for k=1:5]
B=ApproxFun.rrand(Float64,maximum(cols),cols)
cols=Int[rand(1:k+2) for k=1:size(B,1)]
A=ApproxFun.rrand(Float64,maximum(cols),cols)
@test_approx_eq full(A)*full(B) full(A*B)


text:
using KitBase, OffsetArrays, Plots
using KitBase.ProgressMeter: @showprogress

set = KitBase.Setup(
    "scalar", # matter
    "advection", # case
    "1d0f0v", # space
    "gks", # flux
    "", # collision: for scalar conservation laws there are none
    1, # species
    2, # interpolation order
    "vanleer", # limiter
    "period", # boundary
    0.5, # cfl
    1.0, # simulation time
)

pSpace = KitBase.PSpace1D(0.0, 1.0, 100, 1)
vSpace = nothing
property = KitBase.Scalar(0, 1e-6)

w0 = 1.0
prim0 = KitBase.conserve_prim(w0)
ib = nothing
ks = KitBase.SolverSet(set, pSpace, vSpace, property, ib, @__DIR__)

ctr = OffsetArray{KitBase.ControlVolume1D}(undef, eachindex(ks.pSpace.x))
for i in eachindex(ctr)
    u = sin(2π * ks.pSpace.x[i])
    ctr[i] = KitBase.ControlVolume1D(u, KitBase.conserve_prim(u))
end

face = Array{KitBase.Interface1D}(undef, ks.pSpace.nx + 1)
for i = 1:ks.pSpace.nx+1
    face[i] = KitBase.Interface1D(0.0)
end

t = 0.0
dt = KitBase.timestep(ks, ctr, t)
nt = ks.set.maxTime ÷ dt |> Int

anim = @animate for iter = 1:nt
    KitBase.reconstruct!(ks, ctr)

    for i in eachindex(face)
        face[i].fw = KitBase.flux_gks(
            ctr[i-1].w,
            ctr[i].w,
            ks.gas.μᵣ,
            dt,
            0.5 * ks.ps.dx[i-1],
            0.5 * ks.ps.dx[i],
        )
    end

    for i = 1:ks.pSpace.nx
        ctr[i].w += (face[i].fw - face[i+1].fw) / ks.ps.dx[i]
        ctr[i].prim .= KitBase.conserve_prim(ctr[i].w)
    end
    ctr[0].w = ctr[ks.pSpace.nx].w
    ctr[ks.pSpace.nx+1].w = ctr[1].w

    sol = zeros(ks.pSpace.nx)
    for i = 1:ks.pSpace.nx
        sol[i] = ctr[i].w
    end
    plot(ks.pSpace.x[1:ks.pSpace.nx], sol, xlabel = "x", label = "u", ylims = [-1, 1])
end

gif(anim, "burgers.gif", fps = 45)


text:
n = 10
using Test

pairwise_funcs = [
    StatsBase.cov,
    StatsBase.cor,
    StatsBase.countne,
    StatsBase.counteq,
    StatsBase.corkendall,
    StatsBase.corspearman,
    StatsBase.maxad,
    StatsBase.meanad,
    StatsBase.msd,
    StatsBase.psnr,
    StatsBase.rmsd,
    StatsBase.sqL2dist,
    StatsBase.crosscor,
    StatsBase.crosscov
]

@testset "Pairwise statistics on uncertain values" begin
    @testset "$(pairwise_funcs[i])" for i = 1:length(pairwise_funcs)
        f = pairwise_funcs[i]
        @testset for (i, uval) in enumerate(example_uvals)
            if f == StatsBase.psnr
                maxv = 100
                @test f(uval, uval, maxv, n) isa T where T <: Real

            elseif f ∈ [StatsBase.crosscor, StatsBase.crosscov]
                @test f(uval, uval, n) isa AbstractVector{T} where T <: Real
                @test f(uval, uval, 1:5, n) isa AbstractVector{T} where T <: Real
            else
                @test f(uval, uval, n) isa T where T <: Real
            end
        end
    end;
end

# Functions that under the hood use functions with strictly positive domains
# special_pairwise_funcs = [
#     StatsBase.gkldiv,
#     StatsBase.kldivergence,
# ]

# @testset "$(pairwise_funcs[i])" for i = 1:length(pairwise_funcs)
#     f = pairwise_funcs[i]
#     @testset for (i, uval) in enumerate(example_uvals)
#         @test f(uval, uval, n) isa T where T <: Real
#     end
# end;


text:
using DimensionalData, Test

# Tests taken from NamedDims. Thanks @oxinabox

da = DimArray(ones(3), X)
@test Base.BroadcastStyle(typeof(da)) isa DimensionalData.DimensionalStyle

@testset "standard case" begin
    @test da .+ da == 2ones(3)
    @test dims(da .+ da) == dims(da)
    @test da .+ da .+ da == 3ones(3)
    @test dims(da .+ da .+ da) == dims(da)
end

@testset "in place" begin
    @test parent(da .= 1 .* da .+ 7) == 8 * ones(3)
    @test dims(da .= 1 .* da .+ 7) == dims(da)
end

@testset "Dimension disagreement" begin
    @test_throws DimensionMismatch begin
        DimArray(zeros(3, 3, 3), (X, Y, Z)) .+
        DimArray(ones(3, 3, 3), (Y, Z, X))
    end
end

@testset "dims and regular" begin
    da = DimArray(ones(3, 3, 3), (X, Y, Z))
    left_sum = da .+ ones(3, 3, 3)
    @test left_sum == fill(2, 3, 3, 3)
    @test dims(left_sum) == dims(da)
    right_sum = ones(3, 3, 3) .+ da
    @test right_sum == fill(2, 3, 3, 3)
    @test dims(right_sum) == dims(da)
end

@testset "changing type" begin
    @test (da .> 0) isa DimArray
    @test (da .* da .> 0) isa DimArray
    @test (da  .> 0 .> rand(3)) isa DimArray
    @test (da .* rand(3) .> 0.0) isa DimArray
    @test (0 .> da .> 0 .> rand(3)) isa DimArray
    @test (rand(3) .> da  .> 0 .* rand(3)) isa DimArray
    @test (rand(3) .> 1 .> 0 .* da) isa DimArray
end

@testset "broadcasting" begin
    v = DimArray(zeros(3,), X)
    m = DimArray(ones(3, 3), (X, Y))
    s = 0
    @test v .+ m == ones(3, 3) == m .+ v
    @test s .+ m == ones(3, 3) == m .+ s
    @test s .+ v .+ m == ones(3, 3) == m .+ s .+ v
    @test dims(v .+ m) == dims(m .+ v)
    @test dims(s .+ m) == dims(m .+ s)
    @test dims(s .+ v .+ m) == dims(m .+ s .+ v)
end

@testset "adjoint broadcasting" begin
    a = DimArray(reshape(1:12, (4, 3)), (X, Y))
    b = DimArray(1:3, Y)
    @test_throws DimensionMismatch a .* b
    @test_throws DimensionMismatch parent(a) .* parent(b)
    @test parent(a) .* parent(b)' == parent(a .* b')
    @test dims(a .* b') == dims(a)
end

@testset "Mixed array types" begin
    casts = (
        A -> DimArray(A, (X, Y)),  # Named Matrix
        A -> DimArray(A[:, 1], X),  # Named Vector
        A -> DimArray(A[:, 1:1], (X, Y)),  # Named Single Column Matrix
        identity, # Matrix
        A -> A[:, 1], # Vector
        A -> A[:, 1:1], # Single Column Matrix
        first, # Scalar
     )
    for (T1, T2, T3) in Iterators.product(casts, casts, casts)
        all(isequal(identity), (T1, T2, T3)) && continue
        !any(isequal(DimArray), (T1, T2, T3)) && continue
        total = T1(ones(3, 6)) .+ T2(2ones(3, 6)) .+ T3(3ones(3, 6))
        @test total == 6ones(3, 6)
        @test dims(total) == (X(), Y())
    end
end

@testset "in-place assignment .=" begin
    ab = DimArray(rand(2,2), (X, Y))
    ba = DimArray(rand(2,2), (Y, X))
    ac = DimArray(rand(2,2), (X, Z))
    a_ = DimArray(rand(2,2), (X(), AnonDim()))
    z = zeros(2,2)

    @test_throws DimensionMismatch z .= ab .+ ba
    @test_throws DimensionMismatch z .= ab .+ ac
    @test_throws DimensionMismatch a_ .= ab .+ ac
    @test_throws DimensionMismatch ab .= a_ .+ ac
    @test_throws DimensionMismatch ac .= ab .+ ba

    # check that dest is written into:
    @test dims(z .= ab .+ ba') == dims(ab .+ ba')
    @test z == (ab.data .+ ba.data')

    @test dims(z .= ab .+ a_) ==
        (X(Base.OneTo(2); mode=NoIndex()), Y(Base.OneTo(2); mode=NoIndex()))
    @test dims(a_ .= ba' .+ ab) ==
        (X(Base.OneTo(2); mode=NoIndex()), Y(Base.OneTo(2); mode=NoIndex()))
end

@testset "assign using named indexing and dotview" begin
    A = DimArray(zeros(3,2), (X, Y))
    A[X=1:2] .= [1, 2]
    A[X=3] .= 7
    @test A == [1.0 1.0; 2.0 2.0; 7.0 7.0]
end

# @testset "Competing Wrappers" begin
#     da = DimArray(ones(4), X)
#     ta = TrackedArray(5 * ones(4))
#     dt = DimArray(TrackedArray(5 * ones(4)), X)
#     arrays = (da, ta, dt)
#     @testset "$a .- $b"
#     for (a, b) in Iterators.product(arrays, arrays)
#         a === b && continue
#         @test typeof(da .- ta) <: DimArray
#         @test typeof(parent(da .- ta)) <: TrackedArray
#     end
# end


text:
# parallel-random-block-gauss-seidel
function prbgs(A::SparseMatrixCSC, x0::Vector{Float64}, b::Vector{Float64}; block_fraction=0.2, Np=2)

    Nk = length(x0)
    chunksize = Int64(floor(block_fraction*Nk/Np))
    offset = div(chunksize, 2)
    slices = [collect(Iterators.partition(1:Nk, chunksize))...,
              collect(Iterators.partition(offset:Nk, chunksize))...]
    Ns = length(slices)

    x = x0[:]
    r = b - A*x
    bnorm = norm(b)

    err = norm(r) / bnorm
    while err > 1e-10

        picked = [slices[rand(1:end)] for _=1:Np]
        As = []
        for slc in picked
            append!(As, [sparse(view(A, :, slc))])
        end

        deltas = pmap(As) do M
            M \ r
        end

        for i=1:Np
            slc = picked[i]
            x[slc] = x[slc] + deltas[i]
        end

        r[:] = b - A*x
        err = norm(r) / bnorm
    end

    return x
end


text:
function convectionCenteredX(u::AbstractGridGhost, vloc, i, j)
    0.5 * u[i, j] * (u[i+1, j] - u[i-1, j]) / u.grid.dx
        + 0.5 * vloc * (u[i, j+1] - u[i, j-1]) / u.grid.dy
end

function convectionCenteredY(v::AbstractGridGhost, uloc, i, j)
    0.5 * uloc * (v[i+1, j] - v[i-1, j]) / v.grid.dx
        + 0.5 * v[i, j] * (v[i, j+1] - v[i, j-1]) / v.grid.dy
end

function laplacianCentered(g::AbstractGrid, i, j)
    dx2 = 1 / (g.grid.dx * g.grid.dx)
    dy2 = 1 / (g.grid.dy * g.grid.dy)
    c = g[i, j]
    dx2 * (g[i+1, j] - 2. * c + g[i-1, j]) + dy2 * (g[i, j+1] - 2. * c + g[i, j-1])
end

function pressureX(p::AbstractGrid, i, j)
    (p[i, j] - p[i-1, j]) / p.grid.dx
end

function pressureY(p::AbstractGrid, i, j)
    (p[i, j] - p[i, j-1]) / p.grid.dy
end

function fillghost!(g::UGrid, f1, f2)
    ghost1 = g.ghost[1]
    ghost2 = g.ghost[2]
    for i in eachindex(ghost1)
        ghost1[i] = f1(g[i, 1], g[i, 2], g[i, 3])
        ghost2[i] = f2(g[i, end], g[i, end-1], g[i, end-2])
    end
end

function fillghost!(g::VGrid, f1, f2)
    ghost1 = g.ghost[1]
    ghost2 = g.ghost[2]
    for i in eachindex(ghost1)
        ghost1[i] = f1(g[1, i], g[2, i], g[1, i])
        ghost2[i] = f2(g[end, i], g[end-1, i], g[end-2, i])
    end
end

function noslip(arg::Vararg{Real})
    u1, u2, u3 = arg
    -0.2 * (u3 - 5. * u2 + 15. * u1)
end

function naturalbound(arg::Vararg{Real})
    u1 = arg[1]
    u1
end

function veldiv(u::UGrid, v::VGrid, i, j)
    (u[i+1, j] - u[i, j]) / u.grid.dx + (v[i, j+1] - v[i, j]) / v.grid.dy
end

function step_AB2!(sim::Simulation, u_n_1::UGrid, v_n_1::VGrid)
    u_n = copy(sim.u)
    v_n = copy(sim.v)
    nu = getnu(sim.params, sim.u.grid)
    xu = getxs(sim.u)
    yu = getys(sim.u)

    xv = getxs(sim.v)
    yv = getys(sim.v)
    for i in 2:size(sim.u)[1]-1
        for j in 1:size(sim.u)[2]
    # for ci in eachindex(inner_u)
            # (i, j) = Tuple(ci)
            vloc = projectionU(v_n, i, j)
            vloc_n_1 = projectionU(v_n_1, i, j)
            # f1 = - 0.5 * (3. * convectionCenteredX(u_n, vloc, i, j) - convectionCenteredX(u_n_1, vloc_n_1, i, j))
            f1 = - AB2(convectionCenteredX(u_n, vloc, i, j), convectionCenteredX(u_n_1, vloc_n_1, i, j))
            f2 = - pressureX(sim.p, i, j)
            f3 = nu*laplacianCentered(u_n, i, j)

            with_logger(sim.logger) do
                @debug "f1: $f1 - f2: $f2 - f3: $f3"
            end

            mask = inbump(sim.canal, xu[i], yu[j])
            # sim.u[ci] = mask ? (f1 + f2 + f3) / (1 + sim.params.dtdτ) : (f1 + f2 + f3)
            fn1 = u_n[i, j] + sim.dt * (f1 + f2 + f3)
            if mask
                sim.u[i, j] = fn1 / (1 + sim.params.dtdτ * ramp_khi(sim))
            else
                sim.u[i, j] = fn1
            end
            # sim.u[i, j] = mask ? (f1 + f2 + f3) / (1 + sim.params.dtdτ) : (f1 + f2 + f3)
        end
    end

    # for ci in eachindex(inner_v)
    for i in 1:size(sim.v)[1]
        for j in 2:size(sim.v)[2]-1
        # (i, j) = Tuple(ci)
            uloc = projectionV(u_n, i, j)
            uloc_n_1 = projectionV(u_n_1, i, j)
            f1 = - AB2(convectionCenteredY(v_n, uloc, i, j), convectionCenteredY(v_n_1, uloc_n_1, i, j))
            f2 = - pressureY(sim.p, i, j)
            f3 = nu*laplacianCentered(v_n, i, j)

            mask = inbump(sim.canal, xv[i], yv[j])

            fn1 = v_n[i, j] + sim.dt * (f1 + f2 + f3)
            if mask
                sim.v[i, j] = fn1 / (1 + sim.params.dtdτ * ramp_khi(sim))
            else
                sim.v[i, j] = fn1
            end
            # sim.v[ci] = mask ? (f1 + f2 + f3) / (1 + sim.params.dtdτ) : (f1 + f2 + f3)
            # sim.v[i, j] = mask ? (f1 + f2 + f3) / (1 + sim.params.dtdτ) : (f1 + f2 + f3)
        end
    end

    boundaries!(sim)
end

function step_euler!(sim::Simulation)
    u_n = copy(sim.u)
    v_n = copy(sim.v)
    nu = getnu(sim.params, sim.u.grid)
    xu = getxs(sim.u)
    yu = getys(sim.u)

    xv = getxs(sim.v)
    yv = getys(sim.v)
    for i in 2:size(sim.u)[1]-1
        for j in 1:size(sim.u)[2]
    # for ci in eachindex(inner_u)
            # (i, j) = Tuple(ci)
            vloc = projectionU(v_n, i, j)
            f1 = - convectionCenteredX(u_n, vloc, i, j)
            f2 = - pressureX(sim.p, i, j)
            f3 = nu*laplacianCentered(u_n, i, j)

            with_logger(sim.logger) do
                @debug "f1: $f1 - f2: $f2 - f3: $f3"
            end

            mask = inbump(sim.canal, xu[i], yu[j])
            # sim.u[ci] = mask ? (f1 + f2 + f3) / (1 + sim.params.dtdτ) : (f1 + f2 + f3)
            fn1 = u_n[i, j] + sim.dt * (f1 + f2 + f3)
            if mask
                sim.u[i, j] = fn1 / (1 + sim.params.dtdτ * ramp_khi(sim))
            else
                sim.u[i, j] = fn1
            end
            # sim.u[i, j] = mask ? (f1 + f2 + f3) / (1 + sim.params.dtdτ) : (f1 + f2 + f3)
        end
    end

    # for ci in eachindex(inner_v)
    for i in 1:size(sim.v)[1]
        for j in 2:size(sim.v)[2]-1
        # (i, j) = Tuple(ci)
            uloc = projectionV(u_n, i, j)
            f1 = - convectionCenteredY(v_n, uloc, i, j)
            f2 = - pressureY(sim.p, i, j)
            f3 = nu*laplacianCentered(v_n, i, j)

            mask = inbump(sim.canal, xv[i], yv[j])

            fn1 = v_n[i, j] + sim.dt * (f1 + f2 + f3)
            if mask
                sim.v[i, j] = fn1 / (1 + sim.params.dtdτ * ramp_khi(sim))
            else
                sim.v[i, j] = fn1
            end
            # sim.v[ci] = mask ? (f1 + f2 + f3) / (1 + sim.params.dtdτ) : (f1 + f2 + f3)
            # sim.v[i, j] = mask ? (f1 + f2 + f3) / (1 + sim.params.dtdτ) : (f1 + f2 + f3)
        end
    end

    boundaries!(sim)
end

function rightboundary!(sim::Simulation)
    yu = getys(sim.u)
    for col in eachcol(sim.u)
        j = col.indices[2]
        col[end] = col[end] - sim.dt * poiseuille(sim, yu[j]) * (col[end] - col[end-1]) / sim.u.grid.dx
    end
end

function boundaries!(sim::Simulation)
    fillghost!(sim.u, noslip, noslip)
    fillghost!(sim.v, noslip, naturalbound)
    rightboundary!(sim)
end

# TO BE VERIFIED, SEEMS NOT GOOD
function pgauss(p::PGrid, i, j, dx2, dy2)
    right = i != p.grid.m ? p[i+1, j] : -p[i, j]
    left = i != 1 ? p[i-1, j] : p[i, j]
    top = j != p.grid.n ? p[i, j+1] : p[i, j] # DP is 0 at walls
    bot = j != 1 ? p[i, j-1] : p[i, j] # DP is 0 at walls

    (right + left) / dx2 + (top + bot) / dy2
end

function gauss_seidel!(p::PGrid, u::UGrid, v::VGrid, α, dt)
    dx2 = u.grid.dx * u.grid.dx
    dy2 = u.grid.dy * u.grid.dy
    coef = 0.5 * dx2 * dy2 / (dx2 + dy2)
    # innerp = @view p[2:end-1, 2:end-1]
    for ci in eachindex(p)
        i, j = Tuple(ci)
        div = - veldiv(u, v, i, j) / dt
        phit = pgauss(p, i, j, dx2, dy2)
        phistar = coef * ( div + phit )
        p[ci] = α * phistar + (1 - α) * p[ci]
    end
end
gauss_seidel!(sim::Simulation) = gauss_seidel!(sim.p, sim.u, sim.v, sim.params.alpha, sim.dt)

function residual(p::PGrid, u::UGrid, v::VGrid, dt)
    R = 0
    # SHOULD TAKE INTO ACCOUNT BOUNDARIES
    is, js = inrange(p)
    for i in is
        for j in js
            # i, j = Tuple(ci)
            lap = laplacianCentered(p, i, j)
            div = veldiv(u, v, i, j)
            R += (lap - div / dt)^2
        end
    end
    R * p.grid.dx * p.grid.dy
end
function residual(sim::Simulation)
    Rs = residual(sim.p, sim.u, sim.v, sim.dt)
    grid = sim.u.grid
    sim.dt * grid.D / getum(sim.params, grid) * sqrt( 1 / (grid.L * grid.D) * Rs )
end

function step_poisson!(sim::Simulation, tol = 1e-3)
    e = 1000.
    n = 0
    maxn = 1000
    while (abs(e) > tol) && (n < maxn)
        gauss_seidel!(sim)
        e = residual(sim)
        n += 1
    end

    if n == maxn
        with_logger(sim.logger) do
            @warn "Gauss Seidel exceeded $maxn iterations! Aborted"
        end
    end

    with_logger(sim.logger) do
        @info "Poisson equation solved after $n iterations"
    end
end

function poisson_project!(u::UGrid, v::VGrid, p::PGrid, dt)
    ustar = copy(u)
    vstar = copy(v)

    is, js = inrange(u)
    for i in is
        for j in js
            u[i, j] = ustar[i, j] - dt / u.grid.dx * (p[i-1, j] - p[i, j])
        end
    end

    is, js = inrange(v)
    for i in is
        for j in js
            v[i, j] = vstar[i, j] - dt / u.grid.dy * (p[i, j] - p[i, j-1])
        end
    end
end
poisson_project!(sim::Simulation) = poisson_project!(sim.u, sim.v, sim.p, sim.dt)
# function euler!(g::AbstractGridGhost, Δt, f::Function)
#     for ei in eachindex(g)
#         i,j = Tuple(ei)
#         g[ei] = g[ei] + Δt * f(g, i, j)
#     end
# end

function euler(u_n, Δt, f_n)
    u_n + Δt * f_n
end

# function AB2(u_n, Δt, f_n, f_n_1)
#     u_n + 0.5 * Δt * (-f_n_1 + 3. * f_n)
# end

function AB2(f_n, f_n_1)
    0.5 * (-f_n_1 + 3. * f_n)
end

function advectionDiffusionX(g::AbstractGridGhost, i, j)

end

function ramp_khi(ramp::Int, n::Int)
    coef = n / ramp
    coef > 1. ? 1. : coef
end

ramp_khi(sim::Simulation) = ramp_khi(sim.params.ramp, sim.step)


text:
mutable struct BGZFBlockStream
	io::IO # underlying compressed data
	block::Vector{UInt8} # current BGZF block
	#zio # TODO: set type?
	zio::BufferedInputStream{Libz.Source{:inflate,BufferedInputStream{BufferedStreams.EmptyStream}}}

	# TODO: keep track of current virtual offset?
	pos::VirtualOffset
end

function BGZFBlockStream(io::IO, block::Vector{UInt8}=UInt8[])
	block = readblock(io,block)
	BGZFBlockStream(io,block,ZlibInflateInputStream(block),VirtualOffset(0))
end
BGZFBlockStream(filename::AbstractString) = BGZFBlockStream(open(filename))

function readblock(io::IO, block::Vector{UInt8})
	pos = position(io)
	header = zeros(UInt8,12+6)
	readbytes!(io,header,12+6)

	# println("magic: ", header[1:4])
	@assert isequal(header[1:4], UInt8[31,139,8,4])

	bSize = (UInt16(header[12+4+1])) | (UInt16(header[12+4+2])<<8)

	seek(io,pos) # TODO: avoid this seek call
	block = resize!(block,bSize+1)
	readbytes!(io,block,bSize+1)

	# println("bSize+1=", bSize+1)

	block
end

position(s::BGZFBlockStream) = s.pos
eof(s::BGZFBlockStream) = length(s.block) == 28 # assumes the only empty block is at the end

function seek(s::BGZFBlockStream, offs::VirtualOffset)
	# TODO: check if seek is necessary
	seek(s.io, coffset(offs))

	s.block = readblock(s.io, s.block)
	s.zio = ZlibInflateInputStream(s.block)

	# seek in zio (how do I do this in a good way?)
	u = uoffset(offs)
	u > 0 && readbytes(s.zio,u)
	s.pos = offs
	return
end
eofblock(s::BGZFBlockStream) = eof(s.zio)

# assumes eof(s.zio) is false
# TODO: allow reads spanning multiple blocks
function trynextread!(s::BGZFBlockStream,r::BamRead)
	# println("s.pos=", s.pos)
	# println("eofblock(s)=", eofblock(s))
	eofblock(s) && seek(s,s.pos) # initialize next block if needed

	eof(s) && return false

	# nextread!(s.zio,r)
	# s.pos += block_size(r) + 4 # total size of read is block_size + the int32 containing the block_size value
	# # if we reached the end of the block, set virtual offset to start of next block
	# eofblock(s) && (s.pos = VirtualOffset(coffset(s.pos)+length(s.block),0))

	nextread!(s,r)

	# println("after read: s.pos=", s.pos)

	true
end

# TODO: make these support ints that are stored in the intersection between two blocks?
function read(s::BGZFBlockStream,::Type{Int32})
	eofblock(s) && seek(s,s.pos) # initialize next block if needed
	v = read(s.zio,Int32)
	s.pos += 4
	eofblock(s) && (s.pos = VirtualOffset(coffset(s.pos)+length(s.block),0))

	# println(s.pos)

	v
end
function read(s::BGZFBlockStream,::Type{UInt32})
	eofblock(s) && seek(s,s.pos) # initialize next block if needed
	v = read(s.zio,UInt32)
	s.pos += 4
	eofblock(s) && (s.pos = VirtualOffset(coffset(s.pos)+length(s.block),0))

	# println(s.pos)

	v
end

function readbytes!(s::BGZFBlockStream,b::Vector{UInt8},nb)
	eofblock(s) && seek(s,s.pos) # initialize next block if needed

	# println("nb: ",nb)

 	nread = readbytes!(s.zio,b,nb)
 	nb = nb-nread
	s.pos += nread
	eofblock(s) && (s.pos = VirtualOffset(coffset(s.pos)+length(s.block),0))
	# println(s.pos)

	# println("nread: ",nread)
	# println("s.pos: ",s.pos)

 	byteoffs = nread+1
 	totalRead = nread
 	while nb>0
 		eofblock(s) && seek(s,s.pos) # initialize next block if needed

 		# temporary buffer since readbytes! doesn't have API for reading into a buffer with an offset
	 	tempbuf = readbytes(s.zio,nb)
	 	nread = length(tempbuf)

		# println("nb: ",nb)
		# println("nread: ",nread)

	 	b[ byteoffs:byteoffs+nread-1 ] = tempbuf

	 	nb = nb-nread
		s.pos += nread
		eofblock(s) && (s.pos = VirtualOffset(coffset(s.pos)+length(s.block),0))
		# println(s.pos)

	 	byteoffs += nread
	 	totalRead += nread
 	end

 	totalRead
end

# # overwrite read with data from stream
# function nextread!(s::BGZFBlockStream,r::BamRead)
# 	r.block_size = read(io, Int32)
# 	r.refID      = read(io, Int32)
# 	r.pos        = read(io, Int32)
# 	r.bin_mq_nl  = read(io, UInt32)
# 	r.flag_nc    = read(io, UInt32)
# 	r.l_seq      = read(io, Int32)
# 	r.next_refID = read(io, Int32)
# 	r.next_pos   = read(io, Int32)
# 	r.tlen       = read(io, Int32)

# 	reqsize = Int(r.block_size)-8*4
# 	#reqsize > length(r.buf) && error("Read buffer too small!") # we could expand it, but if the file is corrupt, we could allocate huge things here...
# 	@assert reqsize>=0 "Corrupt Bam Read, block_size too small"

# 	# expand read buffer within reasonable limits
# 	# println("pointer: ", pointer(r.buf), ", size: ", length(r.buf))
# 	if reqsize>length(r.buf)
# 		sz = nextpow2(reqsize)
# 		@assert sz<=MAX_READ_SIZE "Read does not fit in maximum read buffer size ($MAX_READ_SIZE bytes)."
# 		r.buf = Array{UInt8,1}(sz)
# 	end

# 	println(reqsize)

# 	n = readbytes!(io,r.buf,reqsize)
# 	n!=reqsize && error("Unexpected end of stream when loading read.")
# 	compute_offsets!(r)
# 	r
# end


text:
include("core.jl")
include("speed.jl")