text:
function layoutable(::Type{<:Colorbar}, fig_or_scene, plot::AbstractPlot; kwargs...)

    layoutable(Colorbar, fig_or_scene;
        colormap = plot.colormap,
        limits = plot.colorrange,
        kwargs...
    )

end

function layoutable(::Type{<:Colorbar}, fig_or_scene, heatmap::Union{Heatmap, Image}; kwargs...)

    layoutable(Colorbar, fig_or_scene;
        colormap = heatmap.colormap,
        limits = heatmap.colorrange,
        highclip = heatmap.highclip,
        lowclip = heatmap.lowclip,
        kwargs...
    )
end

function layoutable(::Type{<:Colorbar}, fig_or_scene, plot::AbstractPlotting.Contourf; kwargs...)

    steps = plot._computed_levels

    limits = lift(steps) do steps
        steps[1], steps[end]
    end

    layoutable(Colorbar, fig_or_scene;
        colormap = plot._computed_colormap,
        limits = limits,
        lowclip = plot._computed_extendlow,
        highclip = plot._computed_extendhigh,
        kwargs...
    )

end

function layoutable(::Type{<:Colorbar}, fig_or_scene; bbox = nothing, kwargs...)
    topscene = get_topscene(fig_or_scene)
    attrs = merge!(Attributes(kwargs), default_attributes(Colorbar, topscene).attributes)

    default_attrs = default_attributes(Colorbar, topscene).attributes
    theme_attrs = subtheme(topscene, :Colorbar)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    @extract attrs (
        label, labelcolor, labelsize, labelvisible, labelpadding, ticklabelsize,
        ticklabelspace, labelfont, ticklabelfont, ticklabelcolor,
        ticklabelsvisible, ticks, tickformat, ticksize, ticksvisible, ticklabelpad, tickalign,
        tickwidth, tickcolor, spinewidth, topspinevisible,
        rightspinevisible, leftspinevisible, bottomspinevisible, topspinecolor,
        leftspinecolor, rightspinecolor, bottomspinecolor, colormap, limits,
        halign, valign, vertical, flipaxisposition, ticklabelalign, flip_vertical_label,
        nsteps, highclip, lowclip,
        minorticksvisible, minortickalign, minorticksize, minortickwidth, minortickcolor, minorticks)

    decorations = Dict{Symbol, Any}()

    protrusions = Node(GridLayoutBase.RectSides{Float32}(0, 0, 0, 0))
    layoutobservables = LayoutObservables{Colorbar}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox, protrusions = protrusions)

    framebox = @lift(round_to_IRect2D($(layoutobservables.computedbbox)))

    highclip_tri_visible = lift(x -> !(isnothing(x) || to_color(x) == to_color(:transparent)), highclip)
    lowclip_tri_visible = lift(x -> !(isnothing(x) || to_color(x) == to_color(:transparent)), lowclip)

    tri_heights = lift(highclip_tri_visible, lowclip_tri_visible, framebox) do hv, lv, box
        if vertical[]
            (lv * width(box), hv * width(box))
        else
            (lv * height(box), hv * height(box))
        end .* sin(pi/3)
    end

    barsize = lift(tri_heights) do heights
        if vertical[]
            max(1, height(framebox[]) - sum(heights))
        else
            max(1, width(framebox[]) - sum(heights))
        end
    end

    barbox = lift(barsize) do sz
        fbox = framebox[]
        if vertical[]
            BBox(left(fbox), right(fbox), bottom(fbox) + tri_heights[][1], top(fbox) - tri_heights[][2])
        else
            BBox(left(fbox) + tri_heights[][1], right(fbox) - tri_heights[][2], bottom(fbox), top(fbox))
        end
    end

    cgradient = lift(colormap, typ = Any) do cmap
        if cmap isa Symbol
            cgrad(cmap)
        elseif cmap isa Tuple{Symbol, Number}
            cgrad(cmap[1], alpha = cmap[2])
        elseif cmap isa PlotUtils.ColorGradient
            cmap
        else
            error("Can't deal with colormap of type $(typeof(cmap))")
        end
    end

    steps = lift(cgradient, nsteps) do cgradient, n
        if cgradient isa PlotUtils.CategoricalColorGradient
            cgradient.values
        else
            collect(LinRange(0, 1, n))
        end::Vector{Float64}
    end

    rects_and_colors = lift(barbox, vertical, steps, cgradient) do bbox, v, steps, gradient
        xmin, ymin = minimum(bbox)
        xmax, ymax = maximum(bbox)

        rects = if v
            yvals = steps .* (ymax - ymin) .+ ymin
            [BBox(xmin, xmax, b, t)
                for (b, t) in zip(yvals[1:end-1], yvals[2:end])]
        else
            xvals = steps .* (xmax - xmin) .+ xmin
            [BBox(l, r, ymin, ymax)
                for (l, r) in zip(xvals[1:end-1], xvals[2:end])]
        end

        colors = get.(Ref(gradient), (steps[1:end-1] .+ steps[2:end]) ./2)

        rects, colors
    end

    rects = poly!(topscene,
        lift(x -> getindex(x, 1), rects_and_colors),
        color = lift(x -> getindex(x, 2), rects_and_colors),
        show_axis = false)

    decorations[:colorrects] = rects

    # hm = heatmap!(topscene, xrange, yrange, colorcells, colormap = colormap, raw = true)
    # decorations[:heatmap] = hm

    highclip_tri = lift(barbox, spinewidth) do box, spinewidth
        if vertical[]
            lb, rb = topline(box)
            l = lb
            r = rb
            t = ((l .+ r) ./ 2) .+ Point2f0(0, sqrt(sum((r .- l) .^ 2)) * sin(pi/3))
            [l, r, t]
        else
            b, t = rightline(box)
            r = ((b .+ t) ./ 2) .+ Point2f0(sqrt(sum((t .- b) .^ 2)) * sin(pi/3), 0)
            [t, b, r]
        end
    end

    highclip_tri_color = Observables.map(highclip) do hc
        to_color(isnothing(hc) ? :transparent : hc)
    end

    highclip_visible = lift(x -> !(isnothing(x) || to_color(x) == to_color(:transparent)), highclip)

    highclip_tri_poly = poly!(topscene, highclip_tri, color = highclip_tri_color,
        strokecolor = :transparent,
        visible = highclip_visible)

    decorations[:highclip] = highclip_tri_poly

    lowclip_tri = lift(barbox, spinewidth) do box, spinewidth
        if vertical[]
            lb, rb = bottomline(box)
            l = lb
            r = rb
            t = ((l .+ r) ./ 2) .- Point2f0(0, sqrt(sum((r .- l) .^ 2)) * sin(pi/3))
            [l, r, t]
        else
            b, t = leftline(box)
            l = ((b .+ t) ./ 2) .- Point2f0(sqrt(sum((t .- b) .^ 2)) * sin(pi/3), 0)
            [b, t, l]
        end
    end

    lowclip_tri_color = Observables.map(lowclip) do lc
        to_color(isnothing(lc) ? :transparent : lc)
    end

    lowclip_visible = lift(x -> !(isnothing(x) || to_color(x) == to_color(:transparent)), lowclip)

    lowclip_tri_poly = poly!(topscene, lowclip_tri, color = lowclip_tri_color,
        strokecolor = :transparent,
        visible = lowclip_visible)

    decorations[:lowclip] = lowclip_tri_poly

    borderpoints = lift(barbox, highclip_visible, lowclip_visible) do bb, hcv, lcv
        if vertical[]
            points = [bottomright(bb), topright(bb)]
            if hcv
                push!(points, highclip_tri[][3])
            end
            append!(points, [topleft(bb), bottomleft(bb)])
            if lcv
                push!(points, lowclip_tri[][3])
            end
            push!(points, bottomright(bb))
            points
        else
            points = [bottomleft(bb), bottomright(bb)]
            if hcv
                push!(points, highclip_tri[][3])
            end
            append!(points, [topright(bb), topleft(bb)])
            if lcv
                push!(points, lowclip_tri[][3])
            end
            push!(points, bottomleft(bb))
            points
        end
    end

    decorations[:spines] = lines!(topscene, borderpoints, linewidth = spinewidth, color = topspinecolor)

    axispoints = lift(barbox, vertical, flipaxisposition) do scenearea,
            vertical, flipaxisposition

        if vertical
            if flipaxisposition
                (bottomright(scenearea), topright(scenearea))
            else
                (bottomleft(scenearea), topleft(scenearea))
            end
        else
            if flipaxisposition
                (topleft(scenearea), topright(scenearea))
            else
                (bottomleft(scenearea), bottomright(scenearea))
            end
        end

    end

    axis = LineAxis(topscene, endpoints = axispoints, flipped = flipaxisposition,
        limits = limits, ticklabelalign = ticklabelalign, label = label,
        labelpadding = labelpadding, labelvisible = labelvisible, labelsize = labelsize,
        labelcolor = labelcolor,
        labelfont = labelfont, ticklabelfont = ticklabelfont, ticks = ticks, tickformat = tickformat,
        ticklabelsize = ticklabelsize, ticklabelsvisible = ticklabelsvisible, ticksize = ticksize,
        ticksvisible = ticksvisible, ticklabelpad = ticklabelpad, tickalign = tickalign,
        tickwidth = tickwidth, tickcolor = tickcolor, spinewidth = spinewidth,
        ticklabelspace = ticklabelspace, ticklabelcolor = ticklabelcolor,
        spinecolor = :transparent, spinevisible = :false, flip_vertical_label = flip_vertical_label,
        minorticksvisible = minorticksvisible, minortickalign = minortickalign,
        minorticksize = minorticksize, minortickwidth = minortickwidth,
        minortickcolor = minortickcolor, minorticks = minorticks)

    decorations[:axis] = axis

    onany(axis.protrusion, vertical, flipaxisposition) do axprotrusion,
            vertical, flipaxisposition

        left, right, top, bottom = 0f0, 0f0, 0f0, 0f0

        if vertical
            if flipaxisposition
                right += axprotrusion
            else
                left += axprotrusion
            end
        else
            if flipaxisposition
                top += axprotrusion
            else
                bottom += axprotrusion
            end
        end

        protrusions[] = GridLayoutBase.RectSides{Float32}(left, right, bottom, top)
    end

    # trigger protrusions with one of the attributes
    vertical[] = vertical[]

    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    Colorbar(fig_or_scene, layoutobservables, attrs, decorations)
end

function tight_ticklabel_spacing!(lc::Colorbar)
    tight_ticklabel_spacing!(lc.elements[:axis])
end


text:
"""
    LaneFollowingAccel
Longitudinal acceleration
"""
immutable LaneFollowingAccel
    a::Float64
end
Base.show(io::IO, a::LaneFollowingAccel) = @printf(io, "LaneFollowingAccel(%6.3f)", a.a)
Base.length(::Type{LaneFollowingAccel}) = 1
Base.convert(::Type{LaneFollowingAccel}, v::Vector{Float64}) = LaneFollowingAccel(v[1])
function Base.copy!(v::Vector{Float64}, a::LaneFollowingAccel)
    v[1] = a.a
    v
end

function propagate(veh::Vehicle1D, action::LaneFollowingAccel, roadway::StraightRoadway, Δt::Float64)

    a = action.a
    s, v = veh.state.s, veh.state.v

    s′ = s + v*Δt + a*Δt*Δt/2
    v′ = v + a*Δt

    s′ = mod_position_to_roadway(s′, roadway)

    return State1D(s′, v′)
end

ext: jl


text:
mydiff(x, y) = mapreduce(abs, +, x - y) / length(x)

run(`unxz -kf sample.raw.xz`)
orig = read_cuboid("sample.raw", 300, 3) |> BitArray

@test mydiff(reconstruct("sample-10.zip"; use_hints = true)[1], orig) < 0.1
@test mydiff(reconstruct("sample-50.zip"; use_hints = true)[1], orig) < 0.05
@test mydiff(reconstruct("sample-90.zip"; use_hints = true)[1], orig) < 0.02


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ErrorEstimationMethod

A method for estimating error of geostatistical solvers.
"""
abstract type ErrorEstimationMethod end

"""
    error(solver, problem, method)

Estimate error of `solver` in a given `problem` with
error estimation `method`.
"""
function error end

# ----------------
# IMPLEMENTATIONS
# ----------------
include("errors/loo.jl")
include("errors/lbo.jl")
include("errors/cv.jl")
include("errors/bcv.jl")
include("errors/wcv.jl")
include("errors/drv.jl")


text:
using Test, MLToolkit

@testset "Nonparametric discrete distribution" begin
    lnpd = LogitNPD()

    k_init = 100
    lnpd1 = LogitNPD(k_init)
    @test length(lnpd1.logitρ) == k_init
    @test lnpd1.logitρ == zeros(k_init)
    @test getlogitρ(lnpd1, 1, k_init) == zeros(k_init)
    @test getρ(lnpd1, 1, k_init) == ones(k_init) / 2
    @test pdf(lnpd1, 0) == 0
    @test pdf(lnpd1, 1) == 0.5
    @test pdf(lnpd1, 2) == 0.25
    @test logpdf(lnpd1, 1) == -log(FT(2))
    @test cdf(lnpd1, 0) == 0
    @test cdf(lnpd1, 1) == 0.5
    @test cdf(lnpd1, 2) == 0.75
    @test invlogcdf(lnpd1, log(FT(0.99))) == 7
    @test ccdf(lnpd1, 0) == 1
    @test ccdf(lnpd1, 1) == 0.5
    @test ccdf(lnpd1, 2) == 0.25
    @test mode(lnpd1) == 1

    k_init = 100
    lnpd1 = LogitNPD(k_init; l_init=1.0)
    @test lnpd1.logitρ == ones(k_init)
    @test getlogitρ(lnpd1, 1, k_init) == ones(k_init)
end


text:
function addmarkers(project;name,i=[],filt=Dict(),colour=:red,markers=:zzz,size=4,label="address")
    # List all existing traces
    tracenames(trace) = trace["name"]
    tracelist = tracenames.(ui.p0.plot.data)
    # Find if called trace exists - if so delete and return, otherwise plot
    i_found = findfirst(tracelist.==name)
    if !isnothing(i_found)
        deletetraces!(ui.p0,i_found)
    else
        if i==[]
            i=1:nrow(project.dat["master"])
        end
        if !isempty(filt)
            findvalue(col,val) = findall(ismissing.(col).==false)[skipmissing(col).==val]
            i = findvalue(project.dat["master"][!,filt["column"]],filt["value"])
        end
        if label=="address"
            labels = project.dat["master"].Full_postal_address[i]
        end

        addtraces!(ui.p0,
            scatter(
                x=project.dat["master"].Easting_HA[i],
                y=project.dat["master"].Northing_HA[i],
                mode="markers",
                marker=attr(color=colour,symbol=markers,size=size),
                line_width=0.1,
                name=name,
                text=labels
            )
        )
    end
end

addmarkers(project;
    name="all",
    colour=:black,
    size=2,
    markers=:O
)

name="before 1919";     addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#1B1464",   size = 7)
name="1919-1929";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#9980FA",   size = 7)
name="1930-1949";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#C4E538",   size = 7)
name="1950-1964";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#EE5A24",   size = 7)
name="1965-1975";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#e1b12c",   size = 7)
name="1976-1983";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#f78fb3",   size = 7)
name="1984-1991";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#3dc1d3",   size = 7)
name="1992-1998";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#5758BB",   size = 7)
name="1999-2002";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#f19066",   size = 7)
name="2003-2007";       addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#c44569",   size = 7)
name="2008 onwards";    addmarkers(project; name=name, filt=Dict("column"=>:EPC_Part_1_Construction_Age_Band,"value"=>name),  colour = "#fed330",   size = 7)

name="before 1919";     addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#1B1464",   size = 7)
name="1919-1929";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#9980FA",   size = 7)
name="1930-1949";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#C4E538",   size = 7)
name="1950-1964";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#EE5A24",   size = 7)
name="1965-1975";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#e1b12c",   size = 7)
name="1976-1983";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#f78fb3",   size = 7)
name="1984-1991";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#3dc1d3",   size = 7)
name="1992-1998";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#5758BB",   size = 7)
name="1999-2002";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#f19066",   size = 7)
name="2003-2007";       addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#c44569",   size = 7)
name="2008 onwards";    addmarkers(project; name=name, filt=Dict("column"=>:ParaDwell_Age_Band,"value"=>name),  colour = "#fed330",   size = 7)


text:
"""
    gpfitpwm(...)

Estimate the GP parameters with the probability weighted moments.

# Implementation

Estimation with the probability weighted moments, as described by [Hosking and Wallis (1987)](https://www.jstor.org/stable/1269343?seq=1),
is only possible in the stationary case.

See also [`gpfitpwm`](@ref) for the other methods, [`gpfit`](@ref), [`gpfitbayes`](@ref) and [`ThresholdExceedance`](@ref).

# Reference

Hosking, J. R. M. and Wallis, J. R. (1987). Parameter and quantile estimation for the Generalized Pareto distribution,
    *Technometrics*, 29:339-349.
"""
function gpfitpwm end

"""
    gpfitpwm(y::Vector{<:Real})

Estimate the GP parameters with the probability weighted moments.

"""
function gpfitpwm(y::Vector{<:Real})::pwmEVA

    model = ThresholdExceedance(Variable("y", y))

    fittedmodel = gpfitpwm(model)

    return fittedmodel

end

"""
    gpfitpwm(df::DataFrame, datacol::Symbol)

Estimate the GP parameters with the probability weighted moments.

Block maxima data are in the column `datacol` of the dataframe `df`.

"""
function gpfitpwm(df::DataFrame, datacol::Symbol)::pwmEVA

    model = ThresholdExceedance(Variable(string(datacol), df[:, datacol]))

    fittedmodel = gpfitpwm(model)

    return fittedmodel

end

"""
    gpfitpwm(model::ThresholdExceedance)

Estimate the GP parameters with the probability weighted moments.

"""
function gpfitpwm(model::ThresholdExceedance)::pwmEVA

    model = validatestationarity(model)

    y = model.data.value

    a₀ = pwm(y,1,0,0)
    a₁ = pwm(y,1,0,1)

    σ̂ = 2*a₀*a₁/(a₀ - 2*a₁)

    k = a₀ / (a₀ - 2*a₁) - 2

    ξ̂ = - k
    ϕ̂ = log(σ̂)

    θ̂ = [ϕ̂; ξ̂]

    fm = pwmEVA{ThresholdExceedance, GeneralizedPareto}(model, θ̂)

    return fm

end


text:
"""
    vcat_axes(x, y) -> Tuple

Returns the appropriate axes for `vcat(x, y)`.

## Examples
```jldoctest
julia> vcat_axes((Index{:a}(1:2), Index{:b}(1:4)), (Index{:z}(1:2), Index(1:4)))
(Index{a}(1:4 => Base.OneTo(4)), Index{b}(1:4 => Base.OneTo(4)))

julia> a, b = [1 2 3 4 5], [6 7 8 9 10; 11 12 13 14 15];

julia> vcat_axes(a, b) == axes(vcat(a, b))
true
```
"""
vcat_axes(x::AbstractArray, y::AbstractArray) = vcat_axes(axes(x), axes(y))
function vcat_axes(x::Tuple, y::Tuple)
    return (cat_axis(first(x), first(y)), combine_indices(tail(x), tail(y))...)
end

"""
    hcat_axes(x, y) -> Tuple

Returns the appropriate axes for `hcat(x, y)`.

## Examples
```jldoctest
julia> hcat_axes((Index{:a}(1:4), Index{:b}(1:2)), (Index{:z}(1:4), Index(1:2)))
(Index{a}(1:4 => Base.OneTo(4)), Index{b}(1:4 => Base.OneTo(4)))

julia> a, b = [1; 2; 3; 4; 5], [6 7; 8 9; 10 11; 12 13; 14 15]

julia> hcat_axes(a, b) == axes(hcat(a, b))
true
```
"""
hcat_axes(x::AbstractArray, y::AbstractArray) = hcat_axes(axes(x), axes(y))
function hcat_axes(x::Tuple, y::Tuple)
    return (combine_index(first(x), first(y)), vcat_axes(tail(x), tail(y))...)
end
hcat_axes(x::NTuple{1,Any}, y::NTuple{1,Any}) = _hcat_axes(combine_index(first(x), first(y)))
_hcat_axes(x) = (x, set_length(unname(x), 2))

"""
    cat_axes(x, y, dims) -> Tuple

Returns the appropriate axes for `cat(x, y; dims)`. If any of `dims` are names
then they should refer to the dimensions of `x`.

## Examples
```jldoctest
julia> cat_axes((Index{:a}(1:4), Index{:b}(1:2)), (Index{:z}(1:4), Index(1:2)), (:a, :b))
(Index{a}(1:8 => Base.OneTo(8)), Index{b}(1:4 => Base.OneTo(4)))
```
"""

"""
    cat_axis(x, y)

Returns the concatenation of the axes `x` and `y`. New subtypes of `AbstractIndex`
must implement a unique `cat_axis` method.
"""
function cat_axis(x::Index, y)
    return Index{cat_names(x, y)}(cat_keys(x, y), cat_values(x, y))
end
function cat_axis(x::SimpleIndex, y)
    return SimpleIndex{cat_names(x, y)}(cat_keys(x, y), cat_values(x, y))
end
cat_axis(x, y) = cat_values(x, y)

"""
    cat_keys(x, y)

Returns the appropriate keys of the `x` and `y` index within the operation `cat_axis(x, y)`

See also: [`cat_axis`](@ref)
"""
cat_keys(x, y) = _cat_keys(keys(x), y)
_cat_keys(x, y) = __cat_keys(Continuity(x), x, y)
__cat_keys(::ContinuousTrait, x, y) = set_length(x, length(x) + length(y))
__cat_keys(::DiscreteTrait, x, y) = make_unique(x, keys(y))

"""
    cat_values(x, y)

Returns the appropriate values of the `x` and `y` index within the operation `cat_axis(x, y)`

See also: [`cat_axis`](@ref)
"""
cat_values(x::AbstractIndex, y) = cat_values(values(x), y)
cat_values(x::AbstractRange, y) = set_length(x, length(x) + length(y))

"""
    cat_names(x, y)

Returns the combined name of `x` and `y` for calls to `cat`. Default behavior is
the same as `combine_names(x, y)`.

See also: [`combine_names`](@ref)
"""
cat_names(x, y) = combine_names(x, y)


text:
using Adapt
using CUDA: CuArray
using OffsetArrays: OffsetArray
using Oceananigans.Fields: fill_halo_regions!
using Oceananigans.Architectures: arch_array

import Oceananigans.TurbulenceClosures: ivd_upper_diagonal,
                                        ivd_lower_diagonal

abstract type AbstractGridFittedBoundary <: AbstractImmersedBoundary end

const GFIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:AbstractGridFittedBoundary}

#####
##### GridFittedBottom (simple 2D immersed boundary)
#####

"""
    GridFittedBottom(bottom)

Return an immersed boundary with an irregular bottom fit to the underlying grid.
"""
struct GridFittedBottom{B} <: AbstractGridFittedBoundary
    bottom :: B
end

const OffsetArrayGridFittedBottom = GridFittedBottom{<:OffsetArray}

function ImmersedBoundaryGrid(grid, ib::OffsetArrayGridFittedBottom)
    TX, TY, TZ = topology(grid)
    return ImmersedBoundaryGrid{TX, TY, TZ}(grid, ib)
end

"""
    ImmersedBoundaryGrid(grid, ib::GridFittedBottom)

Return a grid with `GridFittedBottom` immersed boundary.

Computes ib.bottom and wraps in an array.
"""
function ImmersedBoundaryGrid(grid, ib::GridFittedBottom)
    arch = grid.architecture
    bottom_field = Field{Center, Center, Nothing}(grid)
    set!(bottom_field, ib.bottom)
    fill_halo_regions!(bottom_field)
    offset_bottom_array = dropdims(bottom_field.data, dims=3)
    new_ib = GridFittedBottom(offset_bottom_array)
    return ImmersedBoundaryGrid(grid, new_ib)
end

@inline function is_immersed(i, j, k, underlying_grid, ib::GridFittedBottom)
    x, y, z = node(c, c, c, i, j, k, underlying_grid)
    return @inbounds z < ib.bottom[i, j]
end

on_architecture(arch, ib::GridFittedBottom) = GridFittedBottom(arch_array(arch, ib.bottom))
Adapt.adapt_structure(to, ib::GridFittedBottom) = GridFittedBottom(adapt(to, ib.bottom))

#####
##### Implicit vertical diffusion
#####

####
#### For a center solver we have to check the interface "solidity" at faces k+1 in both the
#### Upper diagonal and the Lower diagonal
#### (because of tridiagonal convention where lower_diagonal on row k is found at k-1)
#### Same goes for the face solver, where we check at centers k in both Upper and lower diagonal
####

@inline immersed_ivd_solid_interface(LX, LY, ::Center, i, j, k, ibg) = solid_interface(LX, LY, Face(), i, j, k+1, ibg)
@inline immersed_ivd_solid_interface(LX, LY, ::Face, i, j, k, ibg)   = solid_interface(LX, LY, Center(), i, j, k, ibg)

# extending the upper and lower diagonal functions of the batched tridiagonal solver

for location in (:upper_, :lower_)
    immersed_func = Symbol(:immersed_ivd_, location, :diagonal)
    func = Symbol(:ivd_ , location, :diagonal)
    @eval begin
        # Disambiguation
        @inline $func(i, j, k, ibg::GFIBG, closure, K, id, LX, LY, LZ::Face, clock, Δt, κz) =
                $immersed_func(i, j, k, ibg::GFIBG, closure, K, id, LX, LY, LZ, clock, Δt, κz)

        @inline $func(i, j, k, ibg::GFIBG, closure, K, id, LX, LY, LZ::Center, clock, Δt, κz) =
                $immersed_func(i, j, k, ibg::GFIBG, closure, K, id, LX, LY, LZ, clock, Δt, κz)

        @inline function $immersed_func(i, j, k, ibg::GFIBG, closure, K, id, LX, LY, LZ, clock, Δt, κz)
            return ifelse(immersed_ivd_solid_interface(LX, LY, LZ, i, j, k, ibg),
                          zero(eltype(ibg.grid)),
                          $func(i, j, k, ibg.grid, closure, K, id, LX, LY, LZ, clock, Δt, κz))
        end
    end
end

#####
##### GridFittedBoundary (experimental 3D immersed boundary)
#####

struct GridFittedBoundary{S} <: AbstractGridFittedBoundary
    mask :: S
end

function ImmersedBoundaryGrid(grid, ib::GridFittedBoundary)
    TX, TY, TZ = topology(grid)
    return ImmersedBoundaryGrid{TX, TY, TZ}(grid, ib)
end

function on_architecture(arch, b::GridFittedBoundary)
    mask = b.mask isa AbstractArray ? arch_array(arch, b.mask) : b.mask
    return GridFittedBoundary(mask)
end

@inline is_immersed(i, j, k, underlying_grid, ib::GridFittedBoundary) = ib.mask(node(c, c, c, i, j, k, underlying_grid)...)


text:
using Juliaup

Juliaup.julia_main()

max_stars_repo_name: piever/juliaup

max_issues_repo_name: piever/juliaup

max_forks_repo_name: piever/juliaup


text:
@testset "helper" begin
    @testset "project function" begin
        @test pso.project(5.0,4.0,6.0,3.0,8.0) == 5.5
        @test pso.project(2.0,0.0,1.0,0.0,1.0) == 2
        @test pso.project(0.3,1.0,0.0,0.0,1.0) == 0.7
    end

    @testset "circle_position function" begin
        @test pso.circle_position(10) ≈ [-1, 0]
        @test pso.circle_position(10, center=[1, 0]) ≈ [0, 0] atol=1E-15
        @test pso.circle_position(10, radius=2, center=[1, 0]) ≈ [-1, 0]
    end

    @testset "swap! function" begin
        vec1 = collect(1:10)
        pso.swap!(vec1, 3, 7)
        @test vec1 == [1, 2, 7, 4, 5, 6, 3, 8, 9, 10]
        vec2 = [1]
        pso.swap!(vec2, 1, 1)
        @test vec2 == [1]
    end
end


text:
#### functions to calculate the imbalance in group representation betwen classes
#### of nodes and the probability of obtaining that imbalance, given a grouping

## Install package if not already installed (slow)
#try Pkg.installed("Iterators") catch Pkg.add("Iterators") end
using Iterators

## Calculate the imbalance of a given class distribution
Imbalance(presvects...) = prod([
        maximum([ vect[ii] for vect in presvects ]) /
        sum([ vect[ii] for vect in presvects ])
    for ii in 1:length(presvects[1]) ])

## Calculate the pvalue of a given imbalance
Pvalue(presvects...) = Float64(prod([
        prod([ binomial(BigInt(sum(vect) - sum(vect[1:(ii - 1)])),
                        BigInt(vect[ii]))
               for vect in presvects ]) /
        binomial(BigInt(sum(sum(presvects)) -
                            sum(sum([ vect[1:(ii - 1)]
                                for vect in presvects ]))),
                 BigInt(sum([ vect[ii] for vect in presvects ])))
    for ii in 1:length(presvects[1]) ]))

## Go through all possible distributions for a given class of species across the
## groups, only keeping those which have the right number of species
function GetPotentialPresences(groupsizes, numspp)
    vects = [0:min(x, numspp) for x in groupsizes]
    presvect = []
    for p in product(vects...)
        if sum(p) != numspp
            continue
        end
        if length(presvect) == 0
            presvect = p
        else
            presvect = [presvect p]
        end
    end
    return presvect
end

## Go through all possible combinations of class distributions, only keeping
## those which do not violate the group sizes (This is the slowest step...)
function GetCompatiblePresences(groupsizes, presvects...)
    ### HARD-CODED ALLOCATION STEP SIZE ###
    STEP = 10000
    ### # # # # # # # # # # # # # # # # ###
    fullpres = Array(Any, STEP)
    counter = 0
    for p in product(presvects...)
        if +([collect(q) for q in p]...) != groupsizes
            continue
        end
#        if length(fullpres) == 0
#            fullpres = [collect(q) for q in p]
#        else
#            fullpres = [fullpres [collect(q) for q in p]]
#        end
        counter += 1
        fullpres[counter] = [collect(q) for q in p]
        if counter == length(fullpres)
            push!(fullpres, Array(Any, STEP))
        end
    end
    return fullpres[1:counter]
end

## Construct a table indicating the pvalue of each possible imbalance value
function GetAllImbalances(presvects)
    imbalances = [Imbalance(presvects[ii]...) for ii in 1:size(presvects, 1)]
    pvalues = [Pvalue(presvects[ii]...) for ii in 1:size(presvects, 1)]
    ## Sanity check
    if sum(pvalues) < .99999999
        error("P-values do not sum to unity: $(round(sum(pvalues), 10))")
    elseif sum(pvalues) > 1.00000001
        error("P-values sum to value greater than unity: $(round(sum(pvalues), 10))")
    end
    return hcat(imbalances, pvalues)
end

## Calculate the pvalue of the empirically observed distribution by summing the
## pvalues of all imbalance values >= to that observed
Emppvalue(empirical, imbatable) = sum(imbatable[Imbalance(empirical...) .<= imbatable[:,1],2])

## Sample several distributions imbalance values to infer the empirical pvalue
function SampleDists(cc, psi, groupvect, classvect, numsamps, outfile)
    imbavals = Vector{Float64}(numsamps)
    uniqClasses = unique(classvect)
    uniqGroups = sort(unique(groupvect))
    randdist = [zeros(Int16, length(uniqGroups)) for xx in 1:length(uniqClasses)]
    for nn in 1:numsamps
        ## equivalent to shuffle groups or classes; groups are integers, so less memory to move
        tmpgroup = shuffle(groupvect)
        for ii in 1:length(uniqClasses)
            for jj in uniqGroups
                ## add one to jj because groups start at 0
                randdist[ii][jj+1] = length(find(yy -> yy == jj,
                                             tmpgroup[find(xx -> xx == uniqClasses[ii],
                                                           classvect)]))
            end
        end
        ## record its imbalance
        imbavals[nn] = Imbalance(randdist...)
        if (numsamps < 100)
            greater = 0
            for imba in imbavals[1:nn]
                if imba >= psi
                    greater = greater + 1
                end
            end
            prob = greater / nn
            write(outfile, "$cc,$psi,$prob,$nn\n")
        elseif (nn % div(numsamps, 100)) == 0
            greater = 0
            for imba in imbavals[1:nn]
                if imba >= psi
                    greater = greater + 1
                end
            end
            prob = greater / nn
            write(outfile, "$cc,$psi,$prob,$nn\n")
        end
    end
    return nothing
end


text:
using GaussianProcesses
using GPR
using Optim
using ConstrainedDynamicsVis
using ConstrainedDynamics
using LineSearches
using Statistics

function _experiment_cp_max(config, id; meandynamics = false)
    traindfs, testdfs = config["datasets"]  # Each thread operates on its own dataset -> no races
    traindf = traindfs.df[id][shuffle(1:nrow(traindfs.df[id]))[1:config["nsamples"]], :]
    testdf = testdfs.df[id][shuffle(1:nrow(testdfs.df[id]))[1:config["testsamples"]], :]
    mechanism = cartpole(1, Δt=0.01, threadlock = config["mechanismlock"])[2]  # Reset Δt to 0.01 in mechanism

    xtest_future_true = [CState(x) for x in testdf.sfuture]
    # Add noise to the dataset
    for df in [traindf, testdf]
        applynoise!(df, config["Σ"], "CP", config["Δtsim"], mechanism.bodies[2].shape.rh[2])
    end
    # Create train and testsets
    xtrain_old = reduce(hcat, [CState(x) for x in traindf.sold])
    xtrain_curr = [CState(x) for x in traindf.scurr]
    vωindices = [9, 22, 23, 24]
    ytrain = [[s[i] for s in xtrain_curr] for i in vωindices]
    xtest_old = [CState(x) for x in testdf.sold]

    predictedstates = Vector{CState{Float64, 2}}()
    params = config["params"]
    gps = Vector{GPE}()
    cache = MDCache()
    for (id, yi) in enumerate(ytrain)
        kernel = SEArd(log.(params[2:end]), log(params[1]))
        mean = meandynamics ? MeanDynamics(mechanism, getμ(vωindices), id, cache) : MeanZero()
        gp = GP(xtrain_old, yi, mean, kernel)
        GaussianProcesses.optimize!(gp, LBFGS(linesearch = BackTracking(order=2)), Optim.Options(time_limit=10.))
        push!(gps, gp)
    end

    projectionerror = 0
    getvω(μ) = return [SVector(0,μ[1],0), SVector(0,μ[2:3]...)], [SVector(zeros(3)...), SVector(μ[4],0,0)]
    for i in 1:length(xtest_old)
        predictedstate, perror = predictdynamics(mechanism, gps, xtest_old[i], config["simsteps"], getvω)
        push!(predictedstates, predictedstate)
        projectionerror += perror
    end
    return predictedstates, xtest_future_true, projectionerror/length(xtest_old)
end

function experiment_cp_mz_max(config, id)
    return _experiment_cp_max(config, id, meandynamics=false)
end

function experiment_cp_md_max(config, id)
    return _experiment_cp_max(config, id, meandynamics=true)
end


text:
# These tests have been adapted from Criterim clojure library

macro fact_max_err(expected, actual, max_err)
    :(@fact abs(($expected) - ($actual)) < $max_err => true)
end

facts("bootstrap estimate") do
    C = Criterion
    testdata = rand(10_000)
    est_mean  = C.bootstrap_bca(testdata, mean, 0.05, 1000)
    @fact (est_mean.point  - 0.5) < 1e-2 => true
    @fact_max_err 0.5 est_mean.lbound  0.02
    @fact_max_err 0.5 est_mean.ubound  0.02
    @fact C.scale_bootstrap_estimate((4.5, 5.5), 0.1) => (0.45, 0.55)
end

facts("bootstrap bca") do
    C = Criterion
    context("mean") do
        test_data = ones(20)
        est_mean  = C.bootstrap_bca(test_data, mean, 0.05, 100)
        @fact est_mean.point  => roughly(1.0)
        @fact est_mean.lbound => roughly(1.0)
        @fact est_mean.ubound => roughly(1.0)
    end

    context("variance") do
        test_data = ones(20)
        est_var   = C.bootstrap_bca(test_data, std, 0.05, 100)
        @fact est_var.point  => roughly(0.0)
        @fact est_var.lbound => roughly(0.0)
        @fact est_var.ubound => roughly(0.0)
    end
end

facts("normal quantile") do
    nq = Criterion.normal_quantile

    @fact nq(0.5001) > 0 => true
    @fact nq(0.4999) < 0 => true
    @fact 2e-8 < (nq(0.999) - nq(0.001)) => true
    let max_error = 1e-7
        @fact_max_err  1.2815515655446   nq(0.9)  max_error
        @fact_max_err  0.674489750196082 nq(0.75) max_error
        @fact_max_err -1.03643338949379  nq(0.15) max_error
        @fact_max_err -2.32634787404084  nq(0.01) max_error
    end
end

facts("normal cdf test") do
    normal_cdf = Criterion.normal_cdf
    max_error = 1.5e-7
    @fact_max_err 0.99865010196837  normal_cdf(3.0) max_error
    @fact_max_err 0.977249868051821 normal_cdf(2.0) max_error
    @fact_max_err 0.841344746068543 normal_cdf(1.0) max_error
    @fact_max_err 0.691462461274013 normal_cdf(0.5) max_error
    @fact_max_err 0.5 normal_cdf(0.0) max_error
    @fact_max_err 0.158655253931457 normal_cdf(-1.0) max_error
    @fact_max_err 0.00134989803163009 normal_cdf(-3.0) max_error
end


text:
function lnn_forward_euler(lnn, lnn_t, a::Assets)

    return a.LAM*(lnn + a.dt*(a.P0*lnn_t + a.P1*a.Sn))
end

function lnTe_forward_euler(lnTe, lnTe_t, a::Assets)

    return a.LAM*(lnTe + a.dt*(a.P0*lnTe_t + a.P1*a.STe))
end

function lnTi_forward_euler(lnTi, lnTi_t, a::Assets)

    return a.LAM*(lnTi + a.dt*(a.P0*lnTi_t + a.P1*a.STi))
end

function u_forward_euler(u, u_t, Te, Ti, a::Assets)

    cs = a.TRGT .* sqrt.(Te + Ti)
    return a.LAM*(u + a.dt*(a.P0*u_t + a.P3*cs))
end

function w_forward_euler(w, w_t, a::Assets)

    return a.LAM*(w + a.dt*a.P0*w_t)
end

function A_forward_euler(A, A_t, Te, Ti, n, phi, er, ev, ad, de2, a::Assets)

    Abohm = a.TRGT .* Abohm_def(Te, Ti, n, phi, er, ev, ad, de2)

    return a.LAM*(A + a.dt*(a.P0*A_t))# + a.P3*Abohm))
end

function diffusion_lnn(lnn, a::Assets)

    rhs = a.P0 * lnn
    return a.DIFF_lnn.A \ rhs
    # return parallel_jacobi_preconditioned_gmres(a.DIFF_lnn, lnn, rhs, 1, a.GRID.Nz)
end

function diffusion_lnTe(lnTe, a::Assets)

    rhs = a.P0 * lnTe
    return a.DIFF_lnTe.A \ rhs
    # return parallel_jacobi_preconditioned_gmres(a.DIFF_lnTe, lnTe, rhs, 1, a.GRID.Nz)
end

function diffusion_lnTi(lnTi, a::Assets)

    rhs = a.P0 * lnTi
    return a.DIFF_lnTi.A \ rhs
    # return parallel_jacobi_preconditioned_gmres(a.DIFF_lnTi, lnTi, rhs, 1, a.GRID.Nz)
end

function diffusion_u(u, Te, Ti, a::Assets)

    cs = a.TRGT .* sqrt.(Te + Ti)
    rhs = a.P0 * u + a.DCHLT3*cs
    return a.DIFF_u.A \ rhs
    # return parallel_jacobi_preconditioned_gmres(a.DIFF_u, u, rhs, 1, a.GRID.Nz)
end

function diffusion_w(w, a::Assets)

    rhs = a.P0 * w
    return a.DIFF_w.A \ rhs
    # return parallel_jacobi_preconditioned_gmres(a.DIFF_w, w, rhs, 1, a.GRID.Nz)
end

function diffusion_A(A, Te, Ti, n, phi, er, ev, ad, de2, a::Assets)

    Abohm = a.TRGT .* Abohm_def(Te, Ti, n, phi, er, ev, ad, de2)
    rhs = a.P0 * A # + a.DCHLT3*Abohm
    return a.DIFF_A.A \ rhs
    # return parallel_jacobi_preconditioned_gmres(a.DIFF_A, A, rhs, 1, a.GRID.Nz)
end

function vorticity_eqn(phi, w, Pi_xx, Pi_yy, Te, n, lnn_x, lnn_y, phi_b, ad, a::Assets)

    M = (a.Dxx + a.Dyy)# + Diagonal(lnn_x)*a.Dx + Diagonal(lnn_y)*a.Dy)
    LHS = LinearLeftHandSide(a.P0*M + a.DCHLT1 + a.DCHLT2 + a.DCHLT3, 2/3)

    lambda = 2.695
    rhs0 = (w - ad*(Pi_xx + Pi_yy)) ./ n
    rhs = a.P0*rhs0 + a.DCHLT1*phi_b + lambda*(a.DCHLT2*Te + a.DCHLT3*Te)
    return LHS.A \ rhs
    # return parallel_jacobi_preconditioned_gmres(LHS, phi, rhs, 1, a.GRID.Nz; m=100)
end

function helmholtz_eqn(A, a::Assets)

    rhs = a.P0 * A
    return a.HHOLTZ.A \ rhs
    # return parallel_jacobi_preconditioned_gmres(a.HHOLTZ, psi, rhs, 1, a.GRID.Nz)
end

function leapfrog!(lnn, lnTe, lnTi, u, w, A, phi, psi, n, Te, Ti, Pe, Pi, j, jn, a::Assets)

    # unpack physical parameters
    am, ad, ki, ke, er, eg, ev, de2, eta = a.params

    # take derivatives of psi
    psi_x = a.Dx*psi
    psi_y = a.Dy*psi
    psi_xy = a.Dxy*psi
    psi_xx = a.Dxx*psi
    psi_yy = a.Dyy*psi

    # subcycle electron thermal conduction terms
    lnTe0 = lnTe[:,2]
    for _=1:a.N_subcycle

        # derivatives
        lnTe_x = a.Dx*lnTe[:,2]
        lnTe_y = a.Dy*lnTe[:,2]
        lnTe_xy = a.Dxy*lnTe[:,2]
        lnTe_xx = a.Dxx*lnTe[:,2]
        lnTe_yy = a.Dyy*lnTe[:,2]
        lnTe_s = partial_s(lnTe[:,2], lnTe_x, lnTe_y, psi_x, psi_y, a.Ds, am)
        lnTe_ss = partial_ss(lnTe[:,2], lnTe_s, lnTe_x, lnTe_y, lnTe_xy, lnTe_xx, lnTe_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, a.Dss, a.Ds, a.Dx, a.Dy, am)

        lnTe_t = lnT_partial_thermal(lnTe_s, lnTe_ss, Te, n, ke)
        lnTe[:,3] = a.LAM *(0.5*(lnTe[:,1] + lnTe[:,2]) + a.dt*(a.P0*lnTe_t + a.P1*a.STe)/a.N_subcycle)

        # derivatives again
        lnTe_x = a.Dx*lnTe[:,2]
        lnTe_y = a.Dy*lnTe[:,2]
        lnTe_xy = a.Dxy*lnTe[:,2]
        lnTe_xx = a.Dxx*lnTe[:,2]
        lnTe_yy = a.Dyy*lnTe[:,2]
        lnTe_s = partial_s(lnTe[:,2], lnTe_x, lnTe_y, psi_x, psi_y, a.Ds, am)
        lnTe_ss = partial_ss(lnTe[:,2], lnTe_s, lnTe_x, lnTe_y, lnTe_xy, lnTe_xx, lnTe_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, a.Dss, a.Ds, a.Dx, a.Dy, am)

        lnTe_t = lnT_partial_thermal(lnTe_s, lnTe_ss, Te, n, ke)
        lnTe[:,3] = a.LAM *(lnTe[:,2] + a.dt*(a.P0*lnTe_t + a.P1*a.STe)/a.N_subcycle)

        # shift indices
        lnTe[:,1:2] = lnTe[:,2:3]
    end
    lnTe[:,2] = 0.5*(lnTe0 + lnTe[:,2])

    # subcycle ion thermal conduction terms
    lnTi0 = lnTi[:,2]
    for _=1:a.N_subcycle

        # derivatives
        lnTi_x = a.Dx*lnTi[:,2]
        lnTi_y = a.Dy*lnTi[:,2]
        lnTi_xy = a.Dxy*lnTi[:,2]
        lnTi_xx = a.Dxx*lnTi[:,2]
        lnTi_yy = a.Dyy*lnTi[:,2]
        lnTi_s = partial_s(lnTi[:,2], lnTi_x, lnTi_y, psi_x, psi_y, a.Ds, am)
        lnTi_ss = partial_ss(lnTi[:,2], lnTi_s, lnTi_x, lnTi_y, lnTi_xy, lnTi_xx, lnTi_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, a.Dss, a.Ds, a.Dx, a.Dy, am)

        lnTi_t = lnT_partial_thermal(lnTi_s, lnTi_ss, Te, n, ke)
        lnTi[:,3] = a.LAM *(0.5*(lnTi[:,1] + lnTi[:,2]) + a.dt*(a.P0*lnTi_t + a.P1*a.STe)/a.N_subcycle)

        # derivatives again
        lnTi_x = a.Dx*lnTi[:,2]
        lnTi_y = a.Dy*lnTi[:,2]
        lnTi_xy = a.Dxy*lnTi[:,2]
        lnTi_xx = a.Dxx*lnTi[:,2]
        lnTi_yy = a.Dyy*lnTi[:,2]
        lnTi_s = partial_s(lnTi[:,2], lnTi_x, lnTi_y, psi_x, psi_y, a.Ds, am)
        lnTi_ss = partial_ss(lnTi[:,2], lnTi_s, lnTi_x, lnTi_y, lnTi_xy, lnTi_xx, lnTi_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, a.Dss, a.Ds, a.Dx, a.Dy, am)

        lnTi_t = lnT_partial_thermal(lnTi_s, lnTi_ss, Te, n, ke)
        lnTi[:,3] = a.LAM *(lnTi[:,2] + a.dt*(a.P0*lnTi_t + a.P1*a.STe)/a.N_subcycle)

        # shift indices
        lnTi[:,1:2] = lnTi[:,2:3]
    end
    lnTi[:,2] = 0.5*(lnTi0 + lnTi[:,2])

    # take other derivatives
    phi_x = a.Dx*phi
    phi_y = a.Dy*phi
    phi_xy = a.Dxy*phi
    phi_xx = a.Dxx*phi
    phi_yy = a.Dyy*phi
    phi_xxx = a.Dxxx*phi
    phi_yyy = a.Dyyy*phi
    phi_xxy = a.Dxxy*phi
    phi_xyy = a.Dxyy*phi
    phi_c = -2*phi_y[:]
    phi_s = partial_s(phi, phi_x, phi_y, psi_x, psi_y, a.Ds, am)

    Pe_x = a.Dx*Pe
    Pe_y = a.Dy*Pe
    Pe_xy = a.Dxy*Pe
    Pe_xx = a.Dxx*Pe
    Pe_yy = a.Dyy*Pe
    Pe_xxx = a.Dxxx*Pe
    Pe_yyy = a.Dyyy*Pe
    Pe_xxy = a.Dxxy*Pe
    Pe_xyy = a.Dxyy*Pe
    Pe_c = -2*Pe_y[:]
    Pe_s = partial_s(Pe, Pe_x, Pe_y, psi_x, psi_y, a.Ds, am)

    Pi_x = a.Dx*Pi
    Pi_y = a.Dy*Pi
    Pi_c = -2*Pi_y[:]
    Pi_s = partial_s(Pi, Pi_x, Pi_y, psi_x, psi_y, a.Ds, am)

    lnn_x = a.Dx*lnn[:,2]
    lnn_y = a.Dy*lnn[:,2]
    lnn_xy = a.Dxy*lnn[:,2]
    lnn_xx = a.Dxx*lnn[:,2]
    lnn_yy = a.Dyy*lnn[:,2]
    lnn_s = partial_s(lnn[:,2], lnn_x, lnn_y, psi_x, psi_y, a.Ds, am)

    lnTe_x = a.Dx*lnTe[:,2]
    lnTe_y = a.Dy*lnTe[:,2]
    lnTe_c = -2*lnTe_y[:]
    lnTe_s = partial_s(lnTe[:,2], lnTe_x, lnTe_y, psi_x, psi_y, a.Ds, am)

    lnTi_x = a.Dx*lnTi[:,2]
    lnTi_y = a.Dy*lnTi[:,2]
    lnTi_c = -2*lnTi_y[:]
    lnTi_s = partial_s(lnTi[:,2], lnTi_x, lnTi_y, psi_x, psi_y, a.Ds, am)

    u_x = a.Dx*u[:,2]
    u_y = a.Dy*u[:,2]
    u_c = -2*u_y[:]
    u_s = partial_s(u[:,2], u_x, u_y, psi_x, psi_y, a.Ds, am)

    j_x = a.Dx*j
    j_y = a.Dy*j
    j_s = partial_s(j, j_x, j_y, psi_x, psi_y, a.Ds, am)

    jn_x = a.Dx*jn
    jn_y = a.Dy*jn
    jn_s = partial_s(jn, jn_x, jn_y, psi_x, psi_y, a.Ds, am)

    G = G_def(Ti, phi_c, Pi_c, u_s, er, ev, ad)
    G_x = a.Dx*G
    G_y = a.Dy*G
    G_c = -2*G_y[:]
    G_s = partial_s(G, G_x, G_y, psi_x, psi_y, a.Ds, am)

    # evaluate total derivatives
    lnn_Dt = lnn_total_ion_derivative(phi_c, Pe_c, n, j_s, u_s, er, ev, ad)
    lnTe_Dt = lnTe_total_electron_derivative_adiabatic(lnn_Dt, j, lnn_s, n, Te, lnTe_c, er, ad)
    lnTi_Dt = lnTi_total_ion_derivative_adiabatic(lnn_Dt, Ti, lnTi_c, er, ad)
    u_Dt = u_total_ion_derivative(Pe_s, Pi_s, G_s, n, Ti, u_c, ev, eg, er, ad)

    # evaluate partial derivatives
    lnn_t = general_partial_derivative(lnn_Dt, lnn_x, lnn_y, lnn_s, phi_x, phi_y, u[:,2])
    lnTe_t = general_partial_derivative(lnTe_Dt, lnTe_x, lnTe_y, lnTe_s, phi_x, phi_y, u[:,2] - jn / ev)
    lnTi_t = general_partial_derivative(lnTi_Dt, lnTi_x, lnTi_y, lnTi_s, phi_x, phi_y, u[:,2])
    u_t = general_partial_derivative(u_Dt, u_x, u_y, u_s, phi_x, phi_y, u[:,2])
    w_t = w_partial_t(n, lnn_x, lnn_y, lnn_xy, lnn_xx, lnn_yy,
                         Pe_x, Pe_y, Pe_xy, Pe_xx, Pe_yy, Pe_xxx, Pe_yyy, Pe_xxy, Pe_xyy, Pe_c,
                         phi_x, phi_y, phi_xy, phi_xx, phi_yy, phi_xxx, phi_yyy, phi_xxy, phi_xyy,
                         G_c, Pi_c, j_s, ad, eg)
    A_t = A_partial_t(phi_x, phi_y, phi_s, j, jn, jn_x, jn_y, jn_s, Pe_s, u[:,2], n, Te, de2, ev, ad, am, eta)

    # compute predictor step
    lnn[:,3] = lnn_forward_euler(0.5*(lnn[:,1] + lnn[:,2]), lnn_t, a)
    lnTe[:,3] = lnTe_forward_euler(0.5*(lnTe[:,1] + lnTe[:,2]), lnTe_t, a)
    lnTi[:,3] = lnTi_forward_euler(0.5*(lnTi[:,1] + lnTi[:,2]), lnTi_t, a)
    u[:,3] = u_forward_euler(0.5*(u[:,1] + u[:,2]), u_t, Te, Ti, a)
    w[:,3] = w_forward_euler(0.5*(w[:,1] + w[:,2]), w_t, a)
    A[:,3] = A_forward_euler(0.5*(A[:,1] + A[:,2]), A_t, Te, Ti, n, phi, er, ev, ad, de2, a)

    # compute inner boundary value for phi
    phi_b = dot(a.FLXAVG, phi) * ones(a.GRID.Nk*a.GRID.Nz)

    # solve linear problems
    lnn[:,3] = diffusion_lnn(lnn[:,3], a)
    lnTe[:,3] = diffusion_lnTe(lnTe[:,3], a)
    lnTi[:,3] = diffusion_lnTi(lnTi[:,3], a)
    u[:,3] = diffusion_u(u[:,3], Te, Ti, a)
    w[:,3] = diffusion_w(w[:,3], a)
    phi[:] = vorticity_eqn(phi, w[:,3], Pi_xx, Pi_yy, Te, n, lnn_x, lnn_y, phi_b, ad, a)
    psi[:] = helmholtz_eqn(A[:,3], a)

    n[:] = exp.(lnn[:,3])
    Te[:] = exp.(lnTe[:,3])
    Ti[:] = exp.(lnTi[:,3])
    Pe[:] = n .* Te
    Pi[:] = n .* Ti
    j[:] = (a.Dxx + a.Dyy) * psi
    jn[:] = j ./ n

    # take derivatives
    psi_x = a.Dx*psi
    psi_y = a.Dy*psi
    psi_xy = a.Dxy*psi
    psi_xx = a.Dxx*psi
    psi_yy = a.Dyy*psi

    phi_x = a.Dx*phi
    phi_y = a.Dy*phi
    phi_xy = a.Dxy*phi
    phi_xx = a.Dxx*phi
    phi_yy = a.Dyy*phi
    phi_xxx = a.Dxxx*phi
    phi_yyy = a.Dyyy*phi
    phi_xxy = a.Dxxy*phi
    phi_xyy = a.Dxyy*phi
    phi_c = -2*phi_y[:]
    phi_s = partial_s(phi, phi_x, phi_y, psi_x, psi_y, a.Ds, am)

    Pe_x = a.Dx*Pe
    Pe_y = a.Dy*Pe
    Pe_xy = a.Dxy*Pe
    Pe_xx = a.Dxx*Pe
    Pe_yy = a.Dyy*Pe
    Pe_xxx = a.Dxxx*Pe
    Pe_yyy = a.Dyyy*Pe
    Pe_xxy = a.Dxxy*Pe
    Pe_xyy = a.Dxyy*Pe
    Pe_c = -2*Pe_y[:]
    Pe_s = partial_s(Pe, Pe_x, Pe_y, psi_x, psi_y, a.Ds, am)

    Pi_x = a.Dx*Pi
    Pi_y = a.Dy*Pi
    Pi_c = -2*Pi_y[:]
    Pi_s = partial_s(Pi, Pi_x, Pi_y, psi_x, psi_y, a.Ds, am)

    lnn_x = a.Dx*lnn[:,3]
    lnn_y = a.Dy*lnn[:,3]
    lnn_xy = a.Dxy*lnn[:,3]
    lnn_xx = a.Dxx*lnn[:,3]
    lnn_yy = a.Dyy*lnn[:,3]
    lnn_s = partial_s(lnn[:,3], lnn_x, lnn_y, psi_x, psi_y, a.Ds, am)

    lnTe_x = a.Dx*lnTe[:,3]
    lnTe_y = a.Dy*lnTe[:,3]
    lnTe_c = -2*lnTe_y[:]
    lnTe_s = partial_s(lnTe[:,3], lnTe_x, lnTe_y, psi_x, psi_y, a.Ds, am)

    lnTi_x = a.Dx*lnTi[:,3]
    lnTi_y = a.Dy*lnTi[:,3]
    lnTi_c = -2*lnTi_y[:]
    lnTi_s = partial_s(lnTi[:,3], lnTi_x, lnTi_y, psi_x, psi_y, a.Ds, am)

    u_x = a.Dx*u[:,3]
    u_y = a.Dy*u[:,3]
    u_c = -2*u_y[:]
    u_s = partial_s(u[:,3], u_x, u_y, psi_x, psi_y, a.Ds, am)

    j_x = a.Dx*j
    j_y = a.Dy*j
    j_s = partial_s(j, j_x, j_y, psi_x, psi_y, a.Ds, am)

    jn_x = a.Dx*jn
    jn_y = a.Dy*jn
    jn_s = partial_s(jn, jn_x, jn_y, psi_x, psi_y, a.Ds, am)

    G = G_def(Ti, phi_c, Pi_c, u_s, er, ev, ad)
    G_x = a.Dx*G
    G_y = a.Dy*G
    G_c = -2*G_y[:]
    G_s = partial_s(G, G_x, G_y, psi_x, psi_y, a.Ds, am)

    # evaluate total derivatives
    lnn_Dt = lnn_total_ion_derivative(phi_c, Pe_c, n, j_s, u_s, er, ev, ad)
    lnTe_Dt = lnTe_total_electron_derivative_adiabatic(lnn_Dt, j, lnn_s, n, Te, lnTe_c, er, ad)
    lnTi_Dt = lnTi_total_ion_derivative_adiabatic(lnn_Dt, Ti, lnTi_c, er, ad)
    u_Dt = u_total_ion_derivative(Pe_s, Pi_s, G_s, n, Ti, u_c, ev, eg, er, ad)

    # evaluate partial derivatives
    lnn_t = general_partial_derivative(lnn_Dt, lnn_x, lnn_y, lnn_s, phi_x, phi_y, u[:,3])
    lnTe_t = general_partial_derivative(lnTe_Dt, lnTe_x, lnTe_y, lnTe_s, phi_x, phi_y, u[:,3] - jn / ev)
    lnTi_t = general_partial_derivative(lnTi_Dt, lnTi_x, lnTi_y, lnTi_s, phi_x, phi_y, u[:,3])
    u_t = general_partial_derivative(u_Dt, u_x, u_y, u_s, phi_x, phi_y, u[:,3])
    w_t = w_partial_t(n, lnn_x, lnn_y, lnn_xy, lnn_xx, lnn_yy,
                         Pe_x, Pe_y, Pe_xy, Pe_xx, Pe_yy, Pe_xxx, Pe_yyy, Pe_xxy, Pe_xyy, Pe_c,
                         phi_x, phi_y, phi_xy, phi_xx, phi_yy, phi_xxx, phi_yyy, phi_xxy, phi_xyy,
                         G_c, Pi_c, j_s, ad, eg)
    A_t = A_partial_t(phi_x, phi_y, phi_s, j, jn, jn_x, jn_y, jn_s, Pe_s, u[:,3], n, Te, de2, ev, ad, am, eta)

    # compute corrector step
    lnn[:,3] = lnn_forward_euler(lnn[:,2], lnn_t, a)
    lnTe[:,3] = lnTe_forward_euler(lnTe[:,2], lnTe_t, a)
    lnTi[:,3] = lnTi_forward_euler(lnTi[:,2], lnTi_t, a)
    u[:,3] = u_forward_euler(u[:,2], u_t, Te, Ti, a)
    w[:,3] = w_forward_euler(w[:,2], w_t, a)
    A[:,3] = A_forward_euler(A[:,2], A_t, Te, Ti, n, phi, er, ev, ad, de2, a)

    # compute inner boundary value for phi
    phi_b = dot(a.FLXAVG, phi) * ones(a.GRID.Nk*a.GRID.Nz)

    # solve linear problems
    lnn[:,3] = diffusion_lnn(lnn[:,3], a)
    lnTe[:,3] = diffusion_lnTe(lnTe[:,3], a)
    lnTi[:,3] = diffusion_lnTi(lnTi[:,3], a)
    u[:,3] = diffusion_u(u[:,3], Te, Ti, a)
    w[:,3] = diffusion_w(w[:,3], a)
    phi[:] = vorticity_eqn(phi, w[:,3], Pi_xx, Pi_yy, Te, n, lnn_x, lnn_y, phi_b, ad, a)
    psi[:] = helmholtz_eqn(A[:,3], a)

    # shift indices
    lnn[:,1:2] = lnn[:,2:3]
    lnTe[:,1:2] = lnTe[:,2:3]
    lnTi[:,1:2] = lnTi[:,2:3]
    u[:,1:2] = u[:,2:3]
    w[:,1:2] = w[:,2:3]
    A[:,1:2] = A[:,2:3]

    n[:] = exp.(lnn[:,3])
    Te[:] = exp.(lnTe[:,3])
    Ti[:] = exp.(lnTi[:,3])
    Pe[:] = n .* Te
    Pi[:] = n .* Ti
    j[:] = (a.Dxx + a.Dyy) * psi
    jn[:] = j ./ n
end


text:
using MeasureInterface
using Test

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


text:
"""
A simple Rewrite Scheduler that applies every rule every time
"""
struct SimpleScheduler <: AbstractScheduler end

cansaturate(s::SimpleScheduler) = true
cansearch(s::SimpleScheduler, r::AbstractRule) = true

function SimpleScheduler(G::EGraph, theory::Vector{<:AbstractRule})
    SimpleScheduler()
end

inform!(s::SimpleScheduler, r, n_matches) = true

setiter!(s::SimpleScheduler, iteration) = nothing


text:
module ApproxInferenceBase
    using Distributions
    using Random
    include("priors.jl")
end # module


text:
# Here, we define additional quantities, gases, units, etc to the
# base one in ThermofluidQuantities

#@displayedunits AreaMomentType      "m^4"    Unitful.𝐋^4
#@displayedunits MomentType          "N*m"    Unitful.𝐌*Unitful.𝐋^2/Unitful.𝐓^2

#@dimvar  Height  LengthType
#@dimvar  Depth   LengthType
@dimvar CentroidX  LengthType
@dimvar CentroidY  LengthType
@dimvar SecondAreaMomentX  AreaMomentType
@dimvar SecondAreaMomentY  AreaMomentType


text:
using Plots
using FFTW

function find_nyquist_frequency(filebase::String, res::Int64)
    arr = Array{Float64}(undef, res)
    for i = 1:res
        for j = 1:res
             arr[j] = sin(j*i*2*pi/res)
        end
        plot(arr, ylim = (-1,1))
        savefig(filebase * string(lpad(string(i-1), 4, string(0))) * ".png")

        plot(abs.(fft(arr)))
        savefig(filebase *"FFT"* string(lpad(string(i-1), 4, string(0)))
                         * ".png")
        println(i)
    end
end

find_nyquist_frequency("images/nyfreq", 40)


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using Behavior.Selection
using Behavior.Selection: TagExpressionInput, SingleTagParser, SequenceParser, Tag

@testset "Selection combinators" begin
    @testset "NotIn" begin
        @testset "Not b; a; OK" begin
            # Arrange
            input = TagExpressionInput("a")
            parser = Selection.NotIn("b")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Char}
            @test result.value == 'a'
        end

        @testset "Not b; c; OK" begin
            # Arrange
            input = TagExpressionInput("c")
            parser = Selection.NotIn("b")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Char}
            @test result.value == 'c'
        end

        @testset "Not b; b; Not OK" begin
            # Arrange
            input = TagExpressionInput("b")
            parser = Selection.NotIn("b")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.BadParseResult{Char}
        end
    end

    @testset "Repeating" begin
        @testset "While not b; b; OK, empty" begin
            # Arrange
            input = TagExpressionInput("b")
            parser = Selection.Repeating(Selection.NotIn("b"))

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Vector{Char}}
            @test result.value == []
        end

        @testset "While not b; ab; OK, a" begin
            # Arrange
            input = TagExpressionInput("ab")
            parser = Selection.Repeating(Selection.NotIn("b"))

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Vector{Char}}
            @test result.value == ['a']
        end

        @testset "While not b; acb; OK, ac" begin
            # Arrange
            input = TagExpressionInput("acb")
            parser = Selection.Repeating(Selection.NotIn("b"))

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Vector{Char}}
            @test result.value == ['a', 'c']
        end
    end

    @testset "TagSelection parser" begin
        @testset "@foo; OK" begin
            # Arrange
            input = TagExpressionInput("@foo")
            parser = SingleTagParser()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.Tag}
            @test result.value == Selection.Tag("@foo")
        end

        @testset "@bar; OK" begin
            # Arrange
            input = TagExpressionInput("@bar")
            parser = SingleTagParser()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.Tag}
            @test result.value == Selection.Tag("@bar")
        end

        @testset "Tag followed by a non-tag character; OK" begin
            nontagchars = "() "

            for nontagchar in nontagchars
                @testset "Non-tag character $(nontagchar)" begin
                    # Arrange
                    input = TagExpressionInput("@foo$(nontagchar)")
                    parser = SingleTagParser()

                    # Act
                    result = parser(input)

                    # Assert
                    @test result isa Selection.OKParseResult{Selection.Tag}
                    @test result.value == Selection.Tag("@foo")
                end
            end
        end
    end

    @testset "TakeUntil parser" begin
        @testset "TakeUntil b; ab; a" begin
            # Arrange
            input = TagExpressionInput("ab")
            parser = Selection.TakeUntil("b")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{String}
            @test result.value == "a"
        end

        @testset "TakeUntil b; aab; aa" begin
            # Arrange
            input = TagExpressionInput("aab")
            parser = Selection.TakeUntil("b")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{String}
            @test result.value == "aa"
        end

        @testset "TakeUntil b; aa; aa" begin
            # Arrange
            input = TagExpressionInput("aa")
            parser = Selection.TakeUntil("b")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{String}
            @test result.value == "aa"
        end

        @testset "TakeUntil space; cc ; cc" begin
            # Arrange
            input = TagExpressionInput("cc ")
            parser = Selection.TakeUntil(" ")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{String}
            @test result.value == "cc"
        end

        @testset "TakeUntil space or ); cc ; cc" begin
            # Arrange
            input = TagExpressionInput("cc ")
            parser = Selection.TakeUntil(" )")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{String}
            @test result.value == "cc"
        end

        @testset "TakeUntil space or ); cc); cc" begin
            # Arrange
            input = TagExpressionInput("cc)")
            parser = Selection.TakeUntil(" )")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{String}
            @test result.value == "cc"
        end

        @testset "TakeUntil b, then TakeUntil d; abcd; a, bc" begin
            # Arrange
            input = TagExpressionInput("abcd")

            # Act
            parser1 = Selection.TakeUntil("b")
            result1 = parser1(input)
            parser2 = Selection.TakeUntil("d")
            result2 = parser2(result1.newinput)

            # Assert
            @test result1 isa Selection.OKParseResult{String}
            @test result1.value == "a"
            @test result2 isa Selection.OKParseResult{String}
            @test result2.value == "bc"
        end

        @testset "TakeUntil space; Prefix whitespace then cc ; cc" begin
            # Arrange
            input = TagExpressionInput(" cc ")
            parser = Selection.TakeUntil(" ")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{String}
            @test result.value == "cc"
        end

        @testset "TakeUntil space twice; a d ; a, then d" begin
            # Arrange
            input = TagExpressionInput("a d ")

            # Act
            parser = Selection.TakeUntil(" ")
            result1 = parser(input)
            result2 = parser(result1.newinput)

            # Assert
            @test result1 isa Selection.OKParseResult{String}
            @test result1.value == "a"
            @test result2 isa Selection.OKParseResult{String}
            @test result2.value == "d"
        end
    end

    @testset "Sequences" begin
        @testset "@foo then @bar; OK" begin
            # Arrange
            input = TagExpressionInput("@foo @bar")
            parser = SequenceParser{Tag}(
                SingleTagParser(),
                SingleTagParser()
            )

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Vector{Tag}}
            @test result.value == [Tag("@foo"), Tag("@bar")]
        end

        @testset "@foo @bar @baz; OK" begin
            # Arrange
            input = TagExpressionInput("@foo @bar @baz")
            parser = SequenceParser{Tag}(
                SingleTagParser(),
                SingleTagParser(),
                SingleTagParser()
            )

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Vector{Tag}}
            @test result.value == [Tag("@foo"), Tag("@bar"), Tag("@baz")]
        end

        @testset "@foo @bar, then a standalone @baz; OK" begin
            # Arrange
            input = TagExpressionInput("@foo @bar @baz")
            parser1 = SequenceParser{Tag}(
                SingleTagParser(),
                SingleTagParser()
            )

            # Act
            result1 = parser1(input)
            parser2 = SingleTagParser()
            result2 = parser2(result1.newinput)

            # Assert
            @test result1 isa Selection.OKParseResult{Vector{Tag}}
            @test result1.value == [Tag("@foo"), Tag("@bar")]
            @test result2 isa Selection.OKParseResult{Tag}
            @test result2.value == Tag("@baz")
        end

        @testset "@a or @b then @c; @c; Not OK" begin
            # Arrange
            input = TagExpressionInput("@c")
            parser = SequenceParser{Selection.TagExpression}(
                Selection.OrParser(),
                SingleTagParser()
            )

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.BadParseResult{Vector{Selection.TagExpression}}
        end
    end

    @testset "Literal parser" begin
        @testset "Literal foo; foo; OK" begin
            # Arrange
            input = TagExpressionInput("foo")
            parser = Selection.Literal("foo")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{String}
            @test result.value == "foo"
        end

        @testset "Literal foo; bar; Not OK" begin
            # Arrange
            input = TagExpressionInput("bar")
            parser = Selection.Literal("foo")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.BadParseResult{String}
        end

        @testset "Literal foobar; bar; Not OK" begin
            # Arrange
            input = TagExpressionInput("bar")
            parser = Selection.Literal("foobar")

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.BadParseResult{String}
        end

        @testset "Literal foo then bar; foobar; OK" begin
            # Arrange
            input = TagExpressionInput("foobar")

            # Act
            parser = SequenceParser{String}(
                Selection.Literal("foo"),
                Selection.Literal("bar")
            )
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Vector{String}}
            @test result.value == ["foo", "bar"]
        end
    end

    @testset "Not parser" begin
        @testset "Not @foo; OK" begin
            # Arrange
            input = TagExpressionInput("not @foo")
            parser = Selection.NotTagParser()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.Not}
            @test result.value == Selection.Not(Selection.Tag("@foo"))
        end
    end

    @testset "Or parser" begin
        @testset "@foo or @bar; OK" begin
            # Arrange
            input = TagExpressionInput("@foo or @bar")
            parser = Selection.OrParser()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.Or}
            @test result.value == Selection.Or(Tag("@foo"), Tag("@bar"))
        end
    end

    @testset "ParenthesesParser" begin
        @testset "Parentheses around tag; (@foo); OK, tag @foo" begin
            # Arrange
            input = TagExpressionInput("(@foo)")
            parser = Selection.ParenthesesParser()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.Parentheses}
            @test result.value == Selection.Parentheses(Tag("@foo"))
        end

        # @testset "Parentheses around Or; (@foo or @bar); OK, @foo or @bar" begin
        #     # Arrange
        #     input = TagExpressionInput("(@foo or @bar)")
        #     parser = Selection.ParenthesesParser()

        #     # Act
        #     result = parser(input)

        #     # Assert
        #     @test result isa Selection.OKParseResult{Selection.Parentheses}
        #     @test result.value == Selection.Parentheses(Selection.Or(Tag("@foo"), Tag("@bar")))
        # end
    end

    @testset "AnyOfParser" begin
        @testset "AnyOf Or,Not; @a or @b; OK, @a or @b" begin
            # Arrange
            input = TagExpressionInput("@a or @b")
            parser = Selection.AnyOfParser(
                Selection.OrParser(),
                Selection.NotTagParser()
            )

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Selection.Or(Tag("@a"), Tag("@b"))
        end

        @testset "AnyOf Or,Not; not @c ; OK, not @c" begin
            # Arrange
            input = TagExpressionInput("not @c")
            parser = Selection.AnyOfParser(
                Selection.OrParser(),
                Selection.NotTagParser()
            )

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Selection.Not(Selection.Tag("@c"))
        end

        @testset "AnyOf Or,Not; @c; Not OK" begin
            # Arrange
            input = TagExpressionInput("@c")
            parser = Selection.AnyOfParser(
                Selection.OrParser(),
                Selection.NotTagParser()
            )

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.BadParseResult{Selection.TagExpression}
        end

        @testset "AnyOf Or,Not,Parentheses; @c; Not OK" begin
            # Arrange
            input = TagExpressionInput("@c")
            parser = Selection.AnyOfParser(
                Selection.OrParser(),
                Selection.NotTagParser(),
                Selection.ParenthesesParser()
            )

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.BadParseResult{Selection.TagExpression}
        end

        @testset "AnyOf Or,Not then @c; not @a @c; OK" begin
            # Arrange
            input = TagExpressionInput("not @a @c")
            parser = Selection.SequenceParser{Selection.TagExpression}(
                Selection.AnyOfParser(
                    Selection.OrParser(),
                    Selection.NotTagParser()
                ),
                SingleTagParser(),
            )

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Vector{Selection.TagExpression}}
            @test result.value == [Selection.Not(Tag("@a")), Tag("@c")]
        end

        @testset "AnyOf Or,Not; @c; AnyOfParser fails, next parser finds @c" begin
            # Arrange
            input = TagExpressionInput("@c")
            parser1 = Selection.AnyOfParser(
                Selection.OrParser(),
                Selection.NotTagParser()
            )
            parser2 = SingleTagParser()

            # Act
            result1 = parser1(input)
            result2 = parser2(result1.newinput)

            # Assert
            @test result2 isa Selection.OKParseResult{Tag}
            @test result2.value == Tag("@c")
        end
    end

    @testset "AnyTagExpression" begin
        @testset "AnyTagExpression; @a; OK" begin
            # Arrange
            input = TagExpressionInput("@a")
            parser = Selection.AnyTagExpression()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Tag("@a")
        end

        @testset "AnyTagExpression; not @a; OK" begin
            # Arrange
            input = TagExpressionInput("not @a")
            parser = Selection.AnyTagExpression()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Selection.Not(Tag("@a"))
        end

        @testset "AnyTagExpression; @a or @b; OK" begin
            # Arrange
            input = TagExpressionInput("@a or @b")
            parser = Selection.AnyTagExpression()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Selection.Or(Tag("@a"), Tag("@b"))
        end

        @testset "AnyTagExpression; (@a); OK" begin
            # Arrange
            input = TagExpressionInput("(@a)")
            parser = Selection.AnyTagExpression()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Selection.Parentheses(Tag("@a"))
        end

        @testset "AnyTagExpression; (@a or @b); OK" begin
            # Arrange
            input = TagExpressionInput("(@a or @b)")
            parser = Selection.AnyTagExpression()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Selection.Parentheses(
                Selection.Or(Tag("@a"), Tag("@b"))
            )
        end

        @testset "AnyTagExpression; (not @c); OK" begin
            # Arrange
            input = TagExpressionInput("(not @c)")
            parser = Selection.AnyTagExpression()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Selection.Parentheses(
                Selection.Not(Tag("@c"))
            )
        end

        @testset "AnyTagExpression; not (@a or @b); OK" begin
            # Arrange
            input = TagExpressionInput("not (@a or @b)")
            parser = Selection.AnyTagExpression()

            # Act
            result = parser(input)

            # Assert
            @test result isa Selection.OKParseResult{Selection.TagExpression}
            @test result.value == Selection.Not(
                Selection.Parentheses(
                    Selection.Or(Tag("@a"), Tag("@b"))
                )
            )
        end
    end
end


text:
using ZChop # For chopping small imaginary parts in ρ
using Distributed

"""
    (t, FI, QFI) = Eff_QFI_HD(Nj, Ntraj, Tfinal, dt; kwargs... )

Evaluate the continuous-time FI and QFI of a final strong measurement for the
estimation of the frequency ω with continuous homodyne monitoring of each half-spin
particle affected by noise at an angle θ, with efficiency η using SME
(stochastic master equation).

The function returns a tuple `(t, FI, QFI)` containing the time vector and the
vectors containing the FI and average QFI

# Arguments

* `Nj`: number of spins
* `Ntraj`: number of trajectories for the SSE
* `Tfinal`: final time of evolution
* `dt`: timestep of the evolution
* `κ = 1`: the noise coupling
* `θ = 0`: noise angle (0 parallel, π/2 transverse)
* `ω = 0`: local value of the frequency
* `η = 1`: measurement efficiency
"""
function Eff_QFI_HD(Nj::Int64,          # Number of spins
    Ntraj::Int64,                       # Number of trajectories
    Tfinal::Float64,                    # Final time
    dt::Float64;                        # Time step
    κ = 1.,                             # Noise coupling
    θ = 0.,                             # Noise angle
    ω = 0.,                             # Frequency
    η = 1.)                             # Measurement efficiency

    Ntime = Int(floor(Tfinal/dt)) # Number of timesteps
    dimJ = Int(2^Nj)   # Dimension of the corresponding Hilbert space

    # Define an array of noise channels
    cj = sqrt(κ/2) *
        [(cos(θ) * σ_j(:z, j, Nj) + sin(θ) * σ_j(:x, j, Nj)) for j = 1:Nj]

    # We store the operators sums for efficiency
    cjSum = [(c + c') for c in cj]

    # We store the operators products for efficiency
    cjProd = [cj[i]*cj[j] for i in eachindex(cj), j in eachindex(cj)]

    # ω is the parameter that we want to estimate
    H = ω * σ(:z, Nj) / 2       # Hamiltonian of the spin system
    dH = σ(:z, Nj) / 2          # Derivative of H wrt the parameter ω

    dW() = sqrt(dt) * randn(Nj) # Define the Wiener increment vector

    # Kraus-like operator, trajectory-independent part

    M0 = sparse(I - 1im * H * dt -
                0.5 * dt * sum([c' * c for c in cj]))

    # Derivative of the Kraus-like operator wrt to ω
    dM = -1im * dH * dt

    # Initial state of the system
    ψ0 = ghz_state(Nj)
    ρ0 = ψ0 * ψ0'

    t = (1 : Ntime) * dt

    # Run evolution for each trajectory, and build up the average
    # for FI and final strong measurement QFI
    result = @distributed (+) for ktraj = 1 : Ntraj
        ρ = ρ0 # Assign initial state to each trajectory

        # Derivative of ρ wrt the parameter
        # Initial state does not depend on the paramter
        dρ = zero(ρ)
        τ = dρ

        # Vectors for the FI and QFI for each trajectory
        FisherT = zero(t)
        QFisherT = zero(t)

        for jt = 1 : Ntime

            # Homodyne current (Eq. 35)
            dy = dt * sqrt(η) * [tr(ρ * c) for c in cjSum] + dW()

            # Kraus operator Eq. (36)
            M = M0 +
                sqrt(η) * sum([cj[j] * dy[j] for j = 1:Nj]) +
                η/2 * sum([cjProd[i,j] *(dy[i] * dy[j] - (i == j ? dt : 0)) for i = 1:Nj, j = 1:Nj])

            # Evolve the density operator
            new_ρ = M * ρ * M' +
                     (1 - η) * dt * sum([c * ρ * c' for c in cj])

            zchop!(new_ρ) # Round off elements smaller than 1e-14

            tr_ρ = tr(new_ρ)

            # Evolve the unnormalized derivative wrt ω
            τ = (M * (τ * M' + ρ * dM') + dM * ρ * M' +
                   (1 - η)* dt * sum([c * τ * c' for c in cj])
                  )/ tr_ρ

            zchop!(τ) # Round off elements smaller than 1e-14

            tr_τ = tr(τ)

            # Now we can renormalize ρ and its derivative wrt ω
            ρ = new_ρ / tr_ρ
            dρ = τ - tr_τ * ρ

            # We evaluate the classical FI for the continuous measurement
            FisherT[jt] = real(tr_τ^2)

            # We evaluate the QFI for a final strong measurement done at time t
            QFisherT[jt] = QFI(ρ, dρ)
        end

        # Use the reduction feature of @distributed for
        # (at the end of each cicle, sum the result to result)
        hcat(FisherT, QFisherT)
    end

    return (t, result[:,1] / Ntraj, result[:,2] / Ntraj)
end


text:
# This file was generated, do not modify it. # hide
r1 = range(xgb, :subsample, lower=0.6, upper=1.0)
r2 = range(xgb, :colsample_bytree, lower=0.6, upper=1.0)
tm = TunedModel(model=xgb, tuning=Grid(resolution=8),
                resampling=CV(rng=234), ranges=[r1,r2],
                measure=cross_entropy)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)


text:
#    Problem : GROUP A
#    *********
#    Extended Cragg and Levy problem.
#    This problem is a sum of m  sets of 5 groups,
#    There are 2m+2 variables. The Hessian matrix is 7-diagonal.
#
#    Original SIF Source:
#    Ph. L. Toint,
#    "Test problems for partially separable optimization and results for the routine
#    PSPMIN", Report 83/4, Department of Mathematics, FUNDP (Namur, B), 1983.
#
#    See also Buckley#18
#
#    Implementation translated from Source:
#    http://eprints.tsu.ge/234/14/Tests%20collection-K-F.pdf
#
#    CRAGGLVY.SIF classification  OUR2-AY-V-0
#
#    Number of variables (variable, based on imput M, the number of group sets)

f = (x) -> begin
    fx = 1.0
    for i in firstindex(x):2:lastindex(x)-3
        item1   = exp(x[i]) - x[i+1]
        item2   = x[i+1] - x[i+2]
        element = x[i+2] - x[i+3]
        item3   = tan(element)+element
        item4   = x[i+3] - 1
        fx     += item1^4 + 100item2^6 + item3^4 + x[i]^8 + item4^2
    end
    return fx
end

g! = (g, x) -> begin
    for i in firstindex(x):2:lastindex(x)-3
        item1   = exp(x[i]) - x[i+1]
        item2   = x[i+1] - x[i+2]
        element = x[i+2] - x[i+3]
        item3   = tan(element)+element
        item4   = x[i+3] - 1
        g[i]   += 4*item1^3*exp(x[i]) + 8*x[i]^7
        g[i+1] += -4item1^3 + 600item2^5
        g[i+2] += -600item2^5 + 4item3^3 * (cos(element)^(-2) + 1)
        g[i+3] += -4.0item3^3*(cos(element)^(-2) + 1.0) + 2*item4
    end
    return g
end

fg! = (g, x) -> begin
    fx = 1.0
    for i in firstindex(x):2:lastindex(x)-3
        item1   = exp(x[i]) - x[i+1]
        item2   = x[i+1] - x[i+2]
        element = x[i+2] - x[i+3]
        item3   = tan(element)+element
        item4   = x[i+3] - 1
        fx     += item1^4 + 100item2^6 + item3^4 + x[i]^8 + item4^2
        g[i]   += 4*item1^3*exp(x[i]) + 8*x[i]^7
        g[i+1] += -4item1^3 + 600item2^5
        g[i+2] += -600item2^5 + 4item3^3 * (cos(element)^(-2) + 1)
        g[i+3] += -4.0item3^3*(cos(element)^(-2) + 1.0) + 2*item4
    end
    return fx, g
end

init = (n::Int=5000) -> begin
	n < 2 && @warn("CRAGGLVY: number of variables must be ≥ 2")
	n = max(n, 2)

    x0 = 2*ones(5000)
    x0[1] = 1.0
    return n, x0
end

@warn "TODO: Debug CRAGGLVY"
# TestSet["CRAGGLVY"] = UncProgram("CRAGGLVY",  f, g!, fg!, init)


text:
module Day15

using AdventOfCode2019
using OffsetArrays

const opposite_direction = Dict(
    0 => 0,
    1 => 2,
    2 => 1,
    3 => 4,
    4 => 3
)

function day15(input::String = readInput(joinpath(@__DIR__, "..", "data", "day15.txt")))
    program = parse.(Int, split(input, ","))
    grid, pos = explore_map(program)
    flood!(grid, pos[1], pos[2], 0)
    return [grid[0, 0], maximum(grid)]
end

function explore_map(program::Array{Int,1})
    inp = Channel{Int}(1)
    out = Channel{Int}(1)
    oxyPos = Channel{CartesianIndex}(1)

    @async AdventOfCode2019.IntCode.run_program!(program, inp, out, nothing)

    gridD = Dict{CartesianIndex{2}, Int}()
    walk!(0, 0, 0, gridD, inp, out, oxyPos)
    pos = take!(oxyPos)

    x1 = minimum(x->x[1] - 1, keys(gridD))
    x2 = maximum(x->x[1] + 1, keys(gridD))
    y1 = minimum(x->x[2] - 1, keys(gridD))
    y2 = maximum(x->x[2] + 1, keys(gridD))
    grid = -1 * ones(Int, x1:x2, y1:y2)
    for (key, entry) in gridD
        grid[key] = entry
    end

    close(inp)
    close(out)
    close(oxyPos)
    return grid, pos
end

function walk!(i::Int, j::Int, backdir::Int, gridD::Dict{CartesianIndex{2}, Int}, inp::Channel{Int}, out::Channel{Int}, oxyPos::Channel{CartesianIndex})
    gridD[CartesianIndex(i, j)] = 0
    for d in 1:4
        di, dj = dir_1d_to_2d(d)
        if !haskey(gridD, CartesianIndex(i + di, j + dj))
            put!(inp, d)
            status = take!(out)
            if status == 0
                gridD[CartesianIndex(i + di, j + dj)] = -2
                continue
            elseif status == 2
                put!(oxyPos, CartesianIndex(i + di, j + dj))
            end
            walk!(i + di, j + dj, opposite_direction[d], gridD, inp, out, oxyPos)
        end
    end
    backdir == 0 && return
    put!(inp, backdir)
    take!(out)
end

function flood!(grid::OffsetArrays.OffsetArray{Int64,2,Array{Int64,2}}, i::Int, j::Int, backdir::Int)
    for d in 1:4
        d == backdir && continue
        di, dj = dir_1d_to_2d(d)
        if grid[i + di, j + dj] > grid[i, j] || grid[i + di, j + dj] == 0
            grid[i + di, j + dj] = grid[i, j] + 1
            flood!(grid, i + di, j + dj, opposite_direction[d])
        end
    end
end

function dir_1d_to_2d(dir::Int)
    dir == 1 && return (-1, 0)
    dir == 2 && return (1, 0)
    dir == 3 && return (0, -1)
    dir == 4 && return (0, 1)
    throw(AssertionError("invalid direction"))
end

end # module


text:
# This file is WIP while the interface for templates is finalized
@testset "Manual Operations Template" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment)
    set_device_model!(template, Line, StaticBranchUnbounded)
    @test !isempty(template.devices)
    @test !isempty(template.branches)
    @test isempty(template.services)
end

@testset "Operations Template Overwrite" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment)
    @test_logs (:info, "Overwriting ThermalStandard existing model") set_device_model!(
        template,
        DeviceModel(ThermalStandard, ThermalBasicUnitCommitment),
    )
    @test template.devices[:ThermalStandard].formulation == ThermalBasicUnitCommitment
end

@testset "Provided Templates Tests" begin
    uc_template = template_unit_commitment()
    @test !isempty(uc_template.devices)
    @test uc_template.devices[:ThermalStandard].formulation == ThermalBasicUnitCommitment
    uc_template = template_unit_commitment(network = DCPPowerModel)
    @test get_transmission_model(uc_template) == DCPPowerModel
    @test !isempty(uc_template.branches)
    @test !isempty(uc_template.services)

    ed_template = template_economic_dispatch()
    @test !isempty(ed_template.devices)
    @test ed_template.devices[:ThermalStandard].formulation == ThermalDispatch
    ed_template = template_economic_dispatch(network = ACPPowerModel)
    @test get_transmission_model(ed_template) == ACPPowerModel
    @test !isempty(ed_template.branches)
    @test !isempty(ed_template.services)
end


text:
module MRIPulse
	#=
		Units philosophy: Use indices in time and frequency domain.
		There are functions to compute the actual units/axes, but use those only
		for verification or plotting.
	=#

	export parallelise, view, normalise!
	export gaussian, sinc
	export index_axis, time_axis, frequency_axis, frequency_unit
	export frequency_shift, frequency_shift!
	export rfspoiling
	export adiabatic_hyperbolic_secant
	export multi_band_pulse, multi_band_pulse!
	export power_integral, peak_power, amplitude_integral
	export minimise_peak_power, grid_evaluate_peak_power # Check with Samy's code (Bioeng KCL)

	import Base: @view, sinc
	import Random
	import Optim
	using LinearAlgebra

	# Macros
	macro view(x, part) # TODO: Extract from here
		if part == :real
			start = :1
		elseif part == :imag
			start = :2
		else
			error("Symbol must be :real or :imag")
		end
		esc(
			quote
				view(
					reinterpret(real(eltype($x)), $x),
					$start:2:2*length($x)
				)
			end
		)
	end

	macro parallelise(threading, forloop) # TODO
		esc(
			quote
				if threading
					@Threads.threads $forloop
				else
					$forloop
				end
			end
		)
	end

	# Pulse shapes
	@inline function gaussian(t::Real, width::Real)::Real # TODO: document what the actual width is
		exp(-(0.5 / width^2) * t^2)
	end

	@inline function sinc(t::Real, width::Real)::Real
		sinc((0.5 / width) * t)
	end

	@inline function rectangular(samples::Integer)::Vector{Float64}
		# Convenience for creating a rectangular pulse that is not breaking the normalise! function
		pulse = Vector{Float64}(undef, samples+2)
		pulse[1] = 0
		pulse[2:end-1] .= 1
		pulse[end] = 0
		return pulse
	end

	# Time and frequency (circular frequency, i.e. in radians) domain Axes
	@inline function index_axis(len::Integer)::UnitRange{Int64}
		@assert len > 0
		(-len ÷ 2):((len + 1) ÷ 2 - 1)
	end

	# Use these for plots only
	@inline function time_axis(len::Integer, δt::R)::StepRangeLen{R} where R <: Real
		@assert δt > 0.0
		δt .* index_axis(len)
	end

	@inline function frequency_unit(len::Integer, δt::Real)::Float64
		@assert len > 0
		@assert δt > 0.0
		2π / (len * δt)
	end

	@inline function frequency_axis(len::Integer, δt::Real)::StepRangeLen{Float64}
		frequency_unit(len, δt) * index_axis(len)
		# Uses twice precision which is weird, follow up
	end

	# Pulse power and flip angle
	@inline function power_integral(pulse::AbstractVector{T})::real(T) where T <: Number
		real(pulse ⋅ pulse)
	end

	function peak_power(pulse::AbstractVector{T})::real(T) where T <: Number
		peak = zero(real(T))
		@inbounds for t ∈ eachindex(pulse)
			peak = max(peak, abs2(pulse[t]))
		end
		return peak
	end

	# Small tip angle approximation
	function amplitude_integral(pulse::AbstractVector{T})::Tuple{real(T), 2} where T <: Number
		amplitude_integral = sum(pulse)
		return abs(amplitude_integral), angle(amplitude_integral)
	end

	# Deapodisation
	@inline function normalise!(x::AbstractArray{T})::AbstractArray{T} where T <: Real
		lower, upper = extrema(x)
		@. x = (x - lower) / (upper - lower)
	end

	# Shift
	@inline function frequency_shift(t::Real, Δω::Real, len::Integer)::ComplexF64
		exp((im * 2π * Δω / len) * t)
	end
	@inline function frequency_shift!(
		pulse::AbstractVector{T} where T <: Number,
		Δω::Real,
		out::V
	)::V where V <: AbstractVector{C} where C <: Complex
		out .= pulse .* frequency_shift.(index_axis(length(pulse)), Δω, length(pulse))
	end
	@inline function frequency_shift(pulse::AbstractArray{T} where T <: Number, Δω::Real)
		frequency_shift!(pulse, Δω, Array{ComplexF64}(undef, length(pulse)))
	end

	# RF spoiling
	function rfspoiling(num::Integer, ϕ0::Real)::Vector{Float64}
		n = 1:num
		ϕ = @. (0.5 * ϕ0) * ( (n - 1) * (n - 2) )
		# Note: At some point it makes sense to compute this differently,
		# since floating point precision decreases.
	end

	# Adiabatic pulses
	macro log_2_sqrt_3()
		log(2 + sqrt(3))
	end
	function adiabatic_hyperbolic_secant(
		t::AbstractVector{<: Real},
		ω1::Real,
		width::Real,
		μ::Real
	)::Tuple{Vector{<: Complex}, Vector{<: Real}}
		#=
			TODO: remove t from args, rather set pulse length T and number samples
			ω1(t) = ω1(0) * sech(βt)
			Δω0(t) = γB0 - Ω1
			where Ω1 is the frequency of the B1 field
			Explain effective field

			Satisfy ellipsoid equation, radii are maximum amplitudes
			1 = ( ω1(t) / ω1(0) )^2 + ( Δω0(t) / Δω0(-∞) )^2

			Solve, setting -μβ = Δω0(t0)
			Δω0(t)	= -μβ * sqrt(1 - ( ω1(t) / ω1(0) )^2 )
					= -μβ * sqrt(1 - sech^2(βt))
					= -μβ * sqrt( (cosh^2(βt) - 1) / cosh^2(βt) )
					= -μβ * sqrt( sinh^2(βt) / cosh^2(βt) )
					= -μβ * tanh(βt) # Silently ignoring the ±

			The phase is then given by
			ϕ(t)	= ∫ Δω0(t') dt' from a given t0 to t
					= -μβ * ∫ tanh(βt') dt' from t0 to t
					= -μ * ∫ tanh(τ) dτ' from βt0 to βt, using τ = βt'
					= -μ * log(cosh(τ)) evaluated at βt0 and βt
					= -μ * ( log(cosh(βt)) - log(cosh(βt0)) )
					= -μ * log(cosh(βt) / cosh(βt0))

			Full width half maximum of the ω1(t) curve

			1/2 = 1 / cosh(βt)
			2 = cosh(βt)
			4 = exp(βt) + exp(-βt)
			exp^2(βt) + 1 - 4 exp(βt) = 0
			exp(βt) = (
				4 ± sqrt(16 - 4 * 1 * 1)
				/ 2
			)
			= 2 ± sqrt(3)

			=> t = 1/β * log(2 ± sqrt(3))
			=> width = 2/β * log(2 + sqrt(3))
			=> β = 2/width * log(2 + sqrt(3))

			Note that lim_{t → ∞} sech(βt) = 2*exp(-βt) can be used to determine a suitable width

			Adiabaticity

			K = |γB_eff / dα/dt|
			= (
				sqrt( ω1^2 * sech^2(βt) + (μβ)^2 * tanh^2(βt))
				/ d/dt arctan(
					ω1 * sech(βt)
					/ (μβ * tanh(βt))
				)
			)
			= (
				sqrt( ω1^2 + (μβ)^2 * sinh^2(βt)) / cosh(βt)
				/ d/dt arctan(
					ω1
					/ (μβ * sinh(βt))
				)
			)
			= (
				sqrt( ω1^2 + (μβ)^2 * sinh^2(βt)) / cosh(βt)
				*
				(1 + (ω1 / μβ)^2 csch^2(βt))
				/
				(ω1 / μ * cosh(βt)/sinh^2(βt))
			)
			= (
				sqrt( ω1^2 + (μβ)^2 * sinh^2(βt))
				*
				(tanh^2(βt) + (ω1 / μβ)^2 sech^2(βt))
				/
				(ω1 / μ)
			)
			= (
				sqrt(1 + (μβ / ω1)^2 * sinh^2(βt))
				*
				(tanh^2(βt) + (ω1 / μβ)^2 sech^2(βt))
				* μ
			)
		=#
		# Compute more convenient width parameter
		β = 2 / width * @log_2_sqrt_3

		# Construct pulse and compute adiabaticity
		pulse = Vector{ComplexF64}(undef, length(t));
		K = Vector{Float64}(undef, length(t));
		let cosh_βt = cosh.(β * t)

			# Pulse
			sech_βt = 1 ./ cosh_βt
			@. pulse = ω1 * sech_βt * exp(-im * μ * log(cosh_βt / cosh(β * t[1])))

			# Adiabaticity
			let η = (μ * β / ω1)^2, tanh_βt = tanh.(β * t)
				@. K = (
					μ * sqrt(1 + η * (cosh_βt^2 - 1))
					* (tanh_βt^2 + 1/η * sech_βt^2)
				)
			end
		end

		return pulse, K
	end

	function adiabatic_phase(ω1::AbstractArray{<: Real}, Δω0_max::Real, dt::Real)
		# Actually, the ellipse equation doesnt have to be fullfilled, I guess it's just easier to check
		# the adibaticity of the effective B-field goes along sth like a circle
		# An adiabatic pulse can be anything, as long as the adibaticity is large enough
		# ω1 must be positive
		Δω0 = Δω0_max .* sqrt.(1 .- (ω1 ./ maximum(ω1)).^2)
		ϕ = cumsum(Δω0) .* dt
		return ϕ
	end

	function adiabaticity(ω1::AbstractVector{<: Real}, Δω0::AbstractVector{<: Real})::Real
		# ω1 = amplitude of B1, Δω0 = offset of B1 rotation to Larmor
		# In units of dt
		tmp = 0
		dα_dt = Vector{Float64}(undef, length(ω1))
		for t = 1:length(ω1)
			dα_dt[t] = arctan(ω1[t] / Δω0[t]) - tmp
			tmp = dα_dt[t]
		end
		A = @. sqrt(ω1^2 + Δω0^2) / abs.(dα_dt)
	end

	# Multi band pulses
	function multi_band_pulse!(
		envelopes::AbstractMatrix{T} where T <: Number,
		Δω::AbstractVector{T} where T <: Real,
		weights::AbstractVector{T} where T <: Number,
		out::V
	)::V where V <: AbstractVector{C} where C <: Complex

		@assert size(envelopes, 1) == length(out)
		@assert size(envelopes, 2) == length(Δω) == length(weights)

		@inbounds for band ∈ eachindex(Δω)
			out .+= (
				weights[band]
				.* view(envelopes, :, band)
				.* frequency_shift.(index_axis(length(out)), Δω[band], length(out))
			)
		end

		return out
	end
	@inline function multi_band_pulse(
		envelopes::AbstractMatrix{T} where T <: Number,
		Δω::AbstractVector{T} where T <: Real,
		weights::AbstractVector{T} where T <: Number
	)
		multi_band_pulse!(envelopes, Δω, weights, zeros(ComplexF64, size(envelopes, 1)))
	end

	# Optimise multi-band pulse phase
	struct OptimisationStorage
		pulse::Vector{ComplexF64}
		expiϕ::Vector{ComplexF64}
	end

	function prepare_optimisation(
		envelopes::AbstractMatrix{T} where T <: Number,
		Δω::AbstractVector{T} where T <: Real,
		weights::AbstractVector{T} where T <: Real,
		threading::Bool
	)::Tuple{Vector{OptimisationStorage}, Vector{ComplexF64}, Matrix{ComplexF64}}

		# Get dimensions
		len, bands = size(envelopes)

		# Compute frequency-shifted pulses, split into the pulse that is invariant
		# and the others for which the phase is modified
		ϕ_const_pulse = weights[1] .* frequency_shift(view(envelopes, :, 1), Δω[1])
		ϕ_var_pulses = zeros(ComplexF64, len, bands - 1)
		@inbounds for band = 2:bands, (i, t_i) ∈ enumerate(index_axis(len))
			ϕ_var_pulses[i, band - 1] = (
				weights[band]
				* envelopes[i, band]
				* frequency_shift(i, Δω[band], len)
			)
		end

		# Pre-allocate memory
		storage = collect(
			OptimisationStorage(zeros(ComplexF64, len), zeros(ComplexF64, bands - 1))
			for _ = 1:(threading ? Threads.nthreads() : 1)
		)

		return storage, ϕ_const_pulse, ϕ_var_pulses
	end

	@inline function objective!(
		ϕ::AbstractVector{Float64},
		storage::OptimisationStorage,
		ϕ_const_pulse::Vector{ComplexF64},
		ϕ_var_pulses::Matrix{ComplexF64}
	)::Float64
		@. storage.expiϕ = exp(im * ϕ)
		mul!(storage.pulse, ϕ_var_pulses, storage.expiϕ)
		@. storage.pulse += ϕ_const_pulse
		peak_power(storage.pulse)
	end

	function minimise_peak_power(
		envelopes::AbstractMatrix{T} where T <: Number,
		Δω::AbstractVector{T} where T <: Real,
		weights::AbstractVector{T} where T <: Real;
		runs::Integer=64,
		threading::Bool=false
	)::Tuple{Vector{ComplexF64}, Vector{Float64}, Float64}

		# Variables used to store the results
		optimised_pulse = nothing
		ϕ = nothing
		minmax_power = Inf

		# Setup pulses, run runs times, store best result in above variables
		let
			( # Type annotations for faster capturing
				storage::Vector{OptimisationStorage},
				ϕ_const_pulse::Vector{ComplexF64},
				ϕ_var_pulses::Matrix{ComplexF64}
			) = prepare_optimisation(envelopes, Δω, weights, threading)

			# rand() will be used to generate initial phases ϕ0;
			# it is not thread-safe need a separate random-number-generator for each
			rng = tuple((Random.MersenneTwister() for _ = 1:Threads.nthreads())...)
			ϕ0 = Vector{Float64}(undef, length(Δω) - 1)

			# Run the optimisation with random initial phases to be more sure about the result
			@parallelise threading for r = 1:runs
				# Initial phase in (0, 2π) elementwise
				@. ϕ0 = π * (
					$rand(rng[Threads.threadid()], 2)
					- $rand(rng[Threads.threadid()], 2) + 1.0
				)
				# Run optimisation
				solution = Optim.optimize(
					ϕ -> objective!(ϕ, storage[Threads.threadid()], ϕ_const_pulse, ϕ_var_pulses),
					0.0, 2π, # Boundaries
					ϕ0,
					Optim.Fminbox(Optim.NelderMead())
				)
				new_minmax_power = Optim.minimum(solution)
				if minmax_power > new_minmax_power
					minmax_power = new_minmax_power
					ϕ = Optim.minimizer(solution)
				end
			end

			# Compute optimised pulse
			let storage = storage[1]
				@. storage.expiϕ = exp(im * ϕ) # TODO: Why is stuff allocated here?
				mul!(storage.pulse, ϕ_var_pulses, storage.expiϕ)
				ϕ_const_pulse .+= storage.pulse
				optimised_pulse = ϕ_const_pulse
			end
		end

		return optimised_pulse, ϕ, minmax_power
	end

	function grid_evaluate_peak_power(
		envelopes::AbstractMatrix{T} where T <: Number,
		Δω::AbstractVector{T} where T <: Real,
		weights::AbstractVector{T} where T <: Real;
		grid::Tuple{Vararg{T, D} where T <: Integer} = (128, 128),
		threading::Bool=false
	)::Array{Float64, D} where D

		max_power = Array{Float64, D}(undef, grid...)
		let
			storage, ϕ_const_pulse, ϕ_var_pulses = (
				prepare_optimisation(envelopes, Δω, weights, threading)
			)

			# Evaluate peak power
			ϕ_axes = collect(range(0, 2π; length=grid[i]) for i = 1:D)
			ϕ = Vector{Float64}(undef, D)
			@parallelise threading for I ∈ CartesianIndices(max_power)
				for d = 1:D
					ϕ[d] = ϕ_axes[d][I[d]]
				end
				max_power[I] = objective!(
					ϕ,
					storage[Threads.threadid()],
					ϕ_const_pulse,
					ϕ_var_pulses
				)
			end
		end

		return max_power
	end

end


text:
"""
This pattern can be apply to subscription behavior.
Suppose we want to subscribe for newspapers.
"""
abstract type Newspaper end
abstract type Subscriber end

function subscribe(::Subscriber, ::Newspaper) end
function unsubscribe(::Subscriber, ::Newspaper) end
function notify(::Newspaper) end

struct SubscriberA <: Subscriber
    name::String
end

update(a::SubscriberA, s::String) = println("$(a.name) is notified by $(s).")

struct SubscriberB <: Subscriber
    name::String
end

update(b::SubscriberB, s::String) = println("$(b.name) is notified by $(s).")

mutable struct AppleNews <: Newspaper
    subscribers::Vector{Subscriber}
end

AppleNews() = AppleNews(Subscriber[])

subscribe(s::Subscriber, a::AppleNews) = push!(a.subscribers, s)
notify(a::AppleNews) = foreach(x -> update(x, "AppleNews"), a.subscribers)

mutable struct BananaNews <: Newspaper
    subscribers::Set{Subscriber}
end

BananaNews() = BananaNews(Set{Subscriber}())

subscribe(s::Subscriber, b::BananaNews) = union!(b.subscribers, [s])
notify(b::BananaNews) = foreach(x -> update(x, "BananaNews"), b.subscribers)