text:
# All functions specific to the Zosterops experiments.
# Initialise predefined species and redefine dispersal and reproduction

# IMPORTANT NOTES (when using `zosterops` mode)
#
# 1. the second environmental niche is now `above-ground carbon` instead of
#    `precipitation` (unfortunately, we can't actually rename it everywhere)
# 2. the `cellsize` setting now determines the patch carrying capacity in
#    individuals, not grams
# 3. the `fertility` setting is the maximum number of offspring per breeding pair per year,
#    instead of a metabolic coefficient. (`noffs = rand(0:fertility)`)
# 4. the `tolerance` setting now determines the probability that a mate of
#    another species is accepted, if no conspecific mate is available
#    (Warning: this only applies when hybridisation == true or speciation == "off")
# 5. `degpleiotropy` must be set to 0, otherwise the species initialisation will fail
# 6. the world must be a rectangle with patch coordinates in row-major order,
#    otherwise dispersal will show weird behaviour (see `coordinate()`)

let zosterops = Individual[] #holds the species archetypes
    """
        initzosteropsspecies()

    Initialise the predefined Zosterops species archetypes.
    """
    function initzosteropsspecies()
        # ensure that "one gene, one trait" is true and that we have species definitions
        (setting("degpleiotropy") != 0) && @simlog("degpleiotropy must be 0", 'e')
        (isnothing(setting("species"))) && @simlog("no species defined", 'e')
        # load per-species trait values from settings (especially AGC optimum and tolerance)
        for sp in setting("species")
            push!(zosterops, initzosteropsspecies(sp))
        end
    end

    """
        initzosteropsspecies(name, precopt, prectol)

    Create a new individual then modify its traits to create a Zosterops archetype.
    """
    function initzosteropsspecies(spectraits::Dict{String,Any})
        archetype = createind(true)
        archetype.lineage = spectraits["lineage"]
        # Find the genes that code for relevant traits and change their values
        for chromosome in archetype.genome
            setting("heterozygosity") && (chromosome.lineage = archetype.lineage)
            for gene in chromosome.genes
                isempty(gene.codes) && continue
                genetrait = setting("traitnames")[gene.codes[1].nameindex]
                if in(genetrait, keys(spectraits))
                    gene.codes[1].value = spectraits[genetrait]
                end
            end
        end
        # then recalculate the individual's trait dict and return the finished archetype
        if setting("linkage") == "none" # `degpleiotropy` is definitely 0 for zosterops mode
            archetype.traits = gettraitdictfast(archetype.genome, setting("traitnames"))
        else
            archetype.traits = gettraitdict(archetype.genome, setting("traitnames"))
        end
        return archetype
    end

    """
        getzosteropsspecies(name, sex)

    Return a new individual of the named species
    """
    global function getzosteropsspecies(name::String, sex::Sex)
        # This function has to be global to escape the let-block of the species list
        isempty(zosterops) && initzosteropsspecies()
        bird = nothing
        for z in zosterops
            if z.lineage == name
                bird = deepcopy(z)
                break
            end
        end
        (isnothing(bird)) && @simlog("Unknown species name: "*name, 'e')
        bird.id = rand(UInt64)
        mutate!(bird, 293.15) #XXX magic number?
        if setting("linkage") == "none" # `degpleiotropy` is definitely 0 for zosterops mode
            bird.traits = gettraitdictfast(bird.genome, setting("traitnames"))
        else
            bird.traits = gettraitdict(bird.genome, setting("traitnames"))
        end
        @simlog("prectol $(bird.traits["prectol"]) precopt $(bird.traits["precopt"]) precadapt $(bird.precadaptation))", 'd')
        bird.size = bird.traits["repsize"] # we need to fix the size after again after mutation
        bird.sex = sex
        return bird
    end
end

"""
    zgenesis(patch)

Create a new community of Zosterops breeding pairs (possibly of multiple species).
Returns an array of individuals.
"""
function zgenesis(patch::Patch)
    community = Array{Individual, 1}()
    (isnothing(setting("species"))) && @simlog("no species defined", 'e')
    # check which species can inhabit this patch
    species = Array{String, 1}()
    for s in setting("species")
        if abs(s["precopt"]-patch.prec) <= s["prectol"]
            push!(species, s["lineage"])
        end
    end
    (isempty(species)) && return community
    # calculate the number of initial breeding pairs and add a male and a female for each
    npairs = Integer(rand(0:round(patch.capacity/2)))
    for i in 1:npairs
        sp = rand(species)
        m = getzosteropsspecies(sp, male)
        f = getzosteropsspecies(sp, female)
        f.partner = m.id
        m.partner = f.id
        push!(community, m)
        push!(community, f)
        @simlog("Adding a pair of Z.$sp", 'd')
    end
    community
end

"""
    zdisperse!(world)

Disperse all juvenile individuals within the world.
"""
function zdisperse!(world::Array{Patch,1})
    for patch in world
        for ind in patch.seedbank
            #growth is one-step, so we integrate it here
            ind.size = ind.traits["repsize"]
            zdisperse!(ind, patch, world)
        end
        patch.seedbank = Array{Individual,1}()
    end
end

"""
    zdisperse!(bird, patch, world)

Dispersal of a single bird. Birds look patches with a suitable
habitat and a free territory or available mate. (Cf. Aben et al. 2016)
"""
function zdisperse!(bird::Individual, patch::Patch, world::Array{Patch,1})
    # keep track of where we've been and calculate the max dispersal distance
    route = [patch.location]
    #XXX disperse!() adds `/sqrt(2)`??
    #XXX sex-biased maximum dispersal?
    !(bird.traits["dispshape"] > 0) && @goto failure # Logistics() requires θ > zero(θ)
    if setting("dispmortality")
        dispmortality= bird.traits["dispmean"]/(setting("dispmean") *setting("dispfactor"))
        rand() < dispmortality && @goto failure
    end
    maxdist = rand(Logistic(bird.traits["dispmean"], bird.traits["dispshape"]))

    while maxdist > 0
        # calculate the best habitat patch in the surroundings (i.e. the closest to AGC optimum)
        bestdest = nothing
        bestfit = nothing
        for pid in shuffle(patch.neighbours) #XXX how expensive is shuffle?
            neighbour = world[pid]
            (neighbour.location in route) && continue
            neighbourfit = abs(neighbour.prec - bird.traits["precopt"])
            if isnothing(bestdest) || neighbourfit < bestfit
                bestdest, bestfit = neighbour, neighbourfit
            end
        end
        (isnothing(bestdest)) && @goto failure
        # check if the patch is within the bird's AGC range and has free space
        if (bestfit <= bird.traits["prectol"] && length(bestdest.community) < bestdest.capacity)
            # can we settle here?
            partner = zfindmate(bestdest.community, bird)
            if !isnothing(partner)
                # if we've found a partner
                bird.partner = partner.id
                partner.partner = bird.id
                @goto success
            elseif length(bestdest.community) <= (bestdest.capacity-2)
                # if there's space for another breeding pair
                @label success
                bird.marked = true
                push!(bestdest.community, bird)
                @simlog("$(idstring(bird)) moved to $(bestdest.location).", 'd')
                return #if we've found a spot, we're done
            end
        end
        push!(route, bestdest.location)
        patch = bestdest
        maxdist -= 1
    end #if the max dispersal distance is reached, the individual simply dies
    @label failure #XXX this could be removed (and `@goto failure` replaced with a simple `return`)
    @simlog("A Z.$(bird.lineage) died after failed dispersal.", 'd')
end

"""
    zfindmate(population, bird)

Find an available breeding partner in the given population, preferring
conspecifics. Returns the partner individual, or nothing.
"""
function zfindmate(population::AbstractArray{Individual, 1}, bird::Individual)
    partner = nothing
    for b in population
        (b.lineage != bird.lineage) && continue
        if ziscompatible(bird, b)
            partner = b
            break
        end
    end
    if isnothing(partner) && (setting("speciation") == "off" || setting("hybridisation") == "on")
        for b in population
            (b.lineage == bird.lineage) && continue
            if ziscompatible(bird, b)
                partner = b
                break
            end
        end
    end
    return partner
end

"""
    ziscompatible(individual1, individual2)

Check to see whether two birds are reproductively compatible.
"""
function ziscompatible(i1::Individual, i2::Individual)
    (i1.sex == i2.sex) && return false
    !(i1.partner == 0 && i2.partner == 0) && return false
    if (setting("speciation") == "off" || setting("hybridisation") == "on") #check for hybridisation
        (i1.lineage != i2.lineage && rand(Float64) > setting("tolerance")) && return false
    else #check for speciation
        (!iscompatible(i1, i2)) && return false
    end
    @simlog("Found a partner: $(idstring(i1)) and $(idstring(i2)).", 'd')
    return true
end

"""
    zreproduce!(patch)

Reproduction of Zosterops breeding pairs in a patch.
"""
function zreproduce!(patch::Patch)
    for bird in patch.community
        @simlog("$(idstring(bird)), $(bird.sex), partner: $(bird.partner)", 'd')
        if bird.partner != 0
            pt = findfirst(b -> b.id == bird.partner, patch.community)
            if isnothing(pt) # birds can find a new partner if one has died
                @simlog("$(idstring(bird)) no longer has a partner.", 'd')
                bird.partner = 0
            elseif bird.sex == female # only mate once per pair
                maxfert = ceil(Int64, setting("fertility")*(bird.precadaptation+0.01))
                noffs = rand(0:maxfert)
                iszero(noffs) && continue
                partner = patch.community[pt]
                @simlog("$(idstring(bird)) mated with $(idstring(partner)).", 'd')
                append!(patch.seedbank, createoffspring(noffs, bird, partner, true))
            end
        end
    end
end

let width = 0, height = 0
    """
        coordinate(x, y, world)

    A utility function to perform a fast look-up for the patch at coordinate x/y.
    Important: this assumes a rectangular world with coordinates in row-major order!
    Returns the index of the desired patch.
    """
    global function coordinate(x::Int, y::Int, world::Array{Patch,1})
        if iszero(width)
            width = maximum(p -> p.location[1], world)
            height = maximum(p -> p.location[2], world)
        end
        (x <= 0 || y <= 0 || x > width || y > height) && return
        i = ((y-1) * width) + x
        return i
    end
end

"""
    findneighbours(world)

Construct a list of neighbours for each patch in the world, for faster lookup
later on. (Must be called during initialisation.)
"""
function findneighbours!(world::Array{Patch,1})
    for patch in world
        for x in (patch.location[1]-1):(patch.location[1]+1)
            for y in (patch.location[2]-1):(patch.location[2]+1)
                ((x,y) == patch.location) && continue
                neighbour = coordinate(x,y,world)
                (isnothing(neighbour)) && continue
                push!(patch.neighbours, neighbour)
            end
        end
    end
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
using Test
using Random
using RandomFourierFeatures
using KernelFunctions
using AbstractGPs

# Writing tests:
# 1. The file structure of the test should match precisely the file structure of src.
#   Amongst other things, this means that there should be exactly 1 test file per src file.
#   This makes it trivially easy for someone to find the tests associated with a particular
#   src file.
# 2. A consequence of 1 is that there should be exactly 1 test file per src file.
# 3. A test file called foo.jl should have the structure:
#   @testset "foo" begin
#       code
#   end
#
#   Note that the testset is called `foo`, not `foo.jl`. Use whatever testset structure
#   seems appropriate within a given file. eg. if multiple types / functions are defined in
#   a particular source file, you might want multiple testsets in the test file.
# 4. Each directory should have its own testset, in which each test file is `include`d.
# 5. Each test file should create its own state, and shouldn't rely on state defined in
#   other test files. If two test files are similar enough to share state, perhaps the
#   corresponding source code should be in the same file.
# 6. If you write a src file without any tests, create a corresponding test file with the
#   usual structure, but without any tests.
# 7. Explicitly create a new random number generate for _at_ _least_ each new test file, and
#   use it whenever generating randomness. This ensures complete control over random number
#   generation and makes it clear what randomness depends on other randomness.
# 8. All `using` statements should appear in runtests.jl.
# 9. List out all test files explicitly (eg. don't loop over them). This makes it easy to
#   disable tests by simply commenting them out, and makes it very clear which tests are not
#   currently being run.
# 10. If utility functionality is required, it should be placed in `src/test_utils.jl` so
#   that other tests can benefit from it.
@info "Packages Loaded"

@testset "RandomFourierFeatures.jl" begin
    include("random_fourier_features.jl")
    println(" ")
    @info "Ran random_fourier_features tests"
end

ext: jl


text:
"""
    Standard()
Standard Scaler. This scaler scale data as:
``\\tilde{\\boldsymbol{x}} = \\frac{x_{i}-\\mu}{\\sigma}``

# Example
```jldoctest preprocessing
julia> x = [
    16.862463771320925 68.10823385851712
    15.382965696961577 65.4313485700859
    8.916228406218375 53.92034559524475
    10.560285659132695 59.17305391117168
    12.142253214135884 62.28708207525656
    5.362107221163482 43.604947901567414
    13.893239446341777 62.44348617377496
    11.871357065173395 60.28433066289655
    29.83792267802442 69.22281924803998
    21.327107214235483 70.15810991597944
    23.852372696012498 69.81780163668844
    26.269031430914108 67.61037566099782
    22.78907104644012 67.78105545358633
    26.73342178134947 68.59263965946904
    9.107259141706415 56.565383817343495
    29.38551885863976 68.1005579469209
    7.935966787763017 53.76264777936664
    29.01677894379809 68.69484161138638
    6.839609488194577 49.69794758177567
    13.95215840314148 62.058116579899085]; #These data are also used to explanations of other functions.

julia> t = [169.80980778351542, 167.9081124078835, 152.30845618985222, 160.3110300206261, 161.96826472170756, 136.02842285615077, 163.98131131382686, 160.117817321485, 172.22758529098235, 172.21342437006865, 171.8939175591617, 169.83018083884602, 171.3878062674257, 170.52487535026015, 156.40282783981309, 170.6488327896672, 151.69267899906185, 172.32478221316322, 145.14365314788827, 163.79383292080666];

julia> scaler = Standard()
Standard(Float64[])

julia> fit!(scaler, x)
2×2 Matrix{Float64}:
 19.0591   64.467
  6.95818   6.68467

julia> transform(scaler, x)
20×2 Matrix{Float64}:
  0.0310374   0.44579
  0.0104337   0.218714
  1.01139     0.710027
  1.37285     0.954091
 -0.895893   -0.270589
  0.983361    0.47397
  1.39855     0.645624
 -1.31861    -0.901482
  0.0147702   0.241708
  0.29675     0.483076
 -0.338824   -0.104812
  0.432442    0.387922
  0.860418    0.567095
 -0.495306    0.140871
  0.963084    0.552767
 -1.38901    -1.05926
  0.469323    0.719196
  0.0669475   0.512023
 -1.47583    -1.44017
 -1.99789    -3.27656
```
"""
mutable struct Standard
    p::AbstractVecOrMat
    Standard() = new(Array{Float64}(undef, 0))
end

"""
    fit!(scaler, x; dims=1)
fit the scaler with `x`. `dims` is the dimension of the number of data.
"""
function fit!(scaler::Standard, x; dims=1)
    scaler.p = vcat(mean(x, dims=dims), std(x, dims=dims))
end

ss(x, m, s) = @. (x-m)/s

"""
    transform(scaler, x; dims=1)
transform data with scaler. `dims` is the dimension of the number of data.
"""
function transform(scaler::Standard, x; dims=1)
    p = scaler.p
    check_size(x, p, dims)
    y = similar(x, Float32)
    if dims == 1
        for i in 1 : size(x, 2)
            y[:, i] = ss(x[:, i], p[:, i]...)
        end
    elseif dims == 2
        for i in 1 : size(x, 1)
            y[i, :] = ss(x[i, :], p[:, i]...)
        end
    end
    return y
end

"""
    fit_transform!(scaler, x; dims=1)
fit scaler with `x`, and transform `x`. `dims` is the dimension of the number of data.
"""
function fit_transform!(scaler::Standard, x; dims=1)
    fit!(scaler, x, dims=dims)
    transform(scaler, x, dims=dims)
end

iss(x, m, s) = @. x*s+m

"""
    inv_transform(scaler, x; dims=1)
Convert `x` in reverse. `dims` is the dimension of the number of data.
"""
function inv_transform(scaler::Standard, x; dims=1)
    p = scaler.p
    check_size(x, p, dims)
    y = similar(x, Float32)
    if dims == 1
        for i in 1 : size(x, 2)
            y[:, i] = iss(x[:, i], p[:, i]...)
        end
    elseif dims == 2
        for i in 1 : size(x, 1)
            y[i, :] = iss(x[i, :], p[:, i]...)
        end
    end
    return y
end


text:
GridLayout(; kwargs...) = GridLayout(1, 1; kwargs...)

observablify(x::Observable) = x
observablify(x, type=Any) = Observable{type}(x)

function GridLayout(nrows::Int, ncols::Int;
        parent = nothing,
        rowsizes = nothing,
        colsizes = nothing,
        addedrowgaps = nothing,
        addedcolgaps = nothing,
        alignmode = Inside(),
        equalprotrusiongaps = (false, false),
        bbox = nothing,
        width::SizeAttribute = Auto(),
        height::SizeAttribute = Auto(),
        tellwidth::Bool = true,
        tellheight::Bool = true,
        halign = :center,
        valign = :center,
        default_rowgap = DEFAULT_ROWGAP_GETTER[](),
        default_colgap = DEFAULT_COLGAP_GETTER[](),
        kwargs...)

    default_rowgap::GapSize = default_rowgap isa Number ? Fixed(default_rowgap)::Fixed : default_rowgap
    default_colgap::GapSize = default_colgap isa Number ? Fixed(default_colgap)::Fixed : default_colgap
    rowsizes = convert_contentsizes(nrows, rowsizes)
    colsizes = convert_contentsizes(ncols, colsizes)
    addedrowgaps = convert_gapsizes(nrows - 1, addedrowgaps, default_rowgap)
    addedcolgaps = convert_gapsizes(ncols - 1, addedcolgaps, default_colgap)

    needs_update = Observable(true)

    content = GridContent[]

    alignmode = observablify(alignmode, AlignMode)
    width = observablify(width)
    height = observablify(height)
    tellwidth = observablify(tellwidth)
    tellheight = observablify(tellheight)
    halign = observablify(halign)
    valign = observablify(valign)

    layoutobservables = layoutobservables = LayoutObservables{GridLayout}(width,
        height, tellwidth, tellheight, halign, valign;
        suggestedbbox = bbox)

    gl = GridLayout(
        parent,
        content, nrows, ncols, rowsizes, colsizes, addedrowgaps,
        addedcolgaps, alignmode, equalprotrusiongaps, needs_update, layoutobservables,
        width, height, tellwidth, tellheight, halign, valign, default_rowgap, default_colgap)

    on(computedbboxobservable(gl)) do cbb
        align_to_bbox!(gl, cbb)
    end

    on(needs_update) do _
        update!(gl)
    end

    gl
end

function update!(gl::GridLayout)
    gl.block_updates && return

    w = determinedirsize(gl, Col())
    h = determinedirsize(gl, Row())

    new_autosize = (w, h)
    new_protrusions = RectSides{Float32}(
        protrusion(gl, Left()),
        protrusion(gl, Right()),
        protrusion(gl, Bottom()),
        protrusion(gl, Top()),
    )

    if autosizeobservable(gl)[] == new_autosize &&
            protrusionsobservable(gl)[] == new_protrusions

        notify(suggestedbboxobservable(gl))
    else
        # otherwise these values will not already be up to date when adding the
        # gridlayout into the next one

        # TODO: this is a double update?
        protrusionsobservable(gl)[] = new_protrusions
        autosizeobservable(gl)[] = new_autosize

        if isnothing(gridcontent(gl))
            notify(suggestedbboxobservable(gl))
        end
    end

    nothing
end

function validategridlayout(gl::GridLayout)
    if gl.nrows < 1
        error("Number of rows can't be smaller than 1")
    end
    if gl.ncols < 1
        error("Number of columns can't be smaller than 1")
    end

    if length(gl.rowsizes) != gl.nrows
        error("There are $nrows rows but $(length(gl.rowsizes)) row sizes.")
    end
    if length(gl.colsizes) != gl.ncols
        error("There are $ncols columns but $(length(gl.colsizes)) column sizes.")
    end
    if length(gl.addedrowgaps) != gl.nrows - 1
        error("There are $nrows rows but $(length(gl.addedrowgaps)) row gaps.")
    end
    if length(gl.addedcolgaps) != gl.ncols - 1
        error("There are $ncols columns but $(length(gl.addedcolgaps)) column gaps.")
    end
end

function with_updates_suspended(f::Function, gl::GridLayout; update = true)
    prev_block_value = gl.block_updates
    gl.block_updates = true
    f()
    gl.block_updates = prev_block_value
    if update
        update!(gl)
    end
    return
end

function connect_layoutobservables!(gc::GridContent)

    disconnect_layoutobservables!(gc::GridContent)

    gc.protrusions_handle = on(protrusionsobservable(gc.content)) do p
        gc.needs_update[] = true
    end
    gc.reportedsize_handle = on(reportedsizeobservable(gc.content)) do c
        gc.needs_update[] = true
    end
end

function disconnect_layoutobservables!(gc::GridContent)
    if !isnothing(gc.protrusions_handle)
        Observables.off(protrusionsobservable(gc.content), gc.protrusions_handle)
        gc.protrusions_handle = nothing
    end
    if !isnothing(gc.reportedsize_handle)
        Observables.off(reportedsizeobservable(gc.content), gc.reportedsize_handle)
        gc.reportedsize_handle = nothing
    end
end

function add_to_gridlayout!(g::GridLayout, gc::GridContent)
    # to be safe
    remove_from_gridlayout!(gc)

    push!(g.content, gc)

    # let the gridcontent know that it's inside a gridlayout
    gc.parent = g
    # change the parent if the gridcontent contains a gridlayout
    content = gc.content
    if content isa GridLayout
        content.parent = g
    end

    on(gc.needs_update) do update
        g.needs_update[] = true
    end

    # trigger relayout
    g.needs_update[] = true
end

function remove_from_gridlayout!(gc::GridContent)
    content = gc.content
    if isnothing(gc.parent)
        if content isa GridLayout
            content.parent = nothing
        end
        return
    end

    i = findfirst(x -> x === gc, gc.parent.content)
    if isnothing(i)
        error("""GridContent had a parent but was not in the content array.
        This must be a bug.""")
    end
    deleteat!(gc.parent.content, i)

    gc.parent = nothing
    # set the parent of a gridlayout content to nothing separately
    # this is mostly for one toplevel parent like a Figure in Makie
    if content isa GridLayout
        content.parent = nothing
    end

    # remove all listeners from needs_update because they could be pointing
    # to previous parents if we're re-nesting layout objects
    empty!(gc.needs_update.listeners)
end

function convert_contentsizes(n, sizes)::Vector{ContentSize}
    if sizes === nothing
        ContentSize[Auto() for _ in 1:n]
    elseif sizes isa ContentSize
        ContentSize[sizes for _ in 1:n]
    elseif sizes isa Vector{<:ContentSize}
        length(sizes) == n ? sizes : error("$(length(sizes)) sizes instead of $n")
    else
        error("Illegal sizes value $sizes")
    end
end

function convert_gapsizes(n, gaps, defaultsize)::Vector{GapSize}
    if gaps === nothing
        GapSize[defaultsize for _ in 1:n]
    elseif gaps isa GapSize
        GapSize[gaps for _ in 1:n]
    elseif gaps isa Vector{<:GapSize}
        length(gaps) == n ? gaps : error("$(length(gaps)) gaps instead of $n")
    else
        error("Illegal gaps value $gaps")
    end
end

function appendrows!(gl::GridLayout, n::Int; rowsizes=nothing, addedrowgaps=nothing, update = true)

    rowsizes = convert_contentsizes(n, rowsizes)
    addedrowgaps = convert_gapsizes(n, addedrowgaps, gl.default_rowgap)

    with_updates_suspended(gl, update = update) do
        gl.nrows += n
        append!(gl.rowsizes, rowsizes)
        append!(gl.addedrowgaps, addedrowgaps)
    end
end

function appendcols!(gl::GridLayout, n::Int; colsizes=nothing, addedcolgaps=nothing, update = true)

    colsizes = convert_contentsizes(n, colsizes)
    addedcolgaps = convert_gapsizes(n, addedcolgaps, gl.default_colgap)

    with_updates_suspended(gl, update = update) do
        gl.ncols += n
        append!(gl.colsizes, colsizes)
        append!(gl.addedcolgaps, addedcolgaps)
    end
end

function prependrows!(gl::GridLayout, n::Int; rowsizes=nothing, addedrowgaps=nothing, update = true)

    rowsizes = convert_contentsizes(n, rowsizes)
    addedrowgaps = convert_gapsizes(n, addedrowgaps, gl.default_rowgap)

    foreach(gl.content) do gc
        span = gc.span
        newspan = Span(span.rows .+ n, span.cols)
        gc.span = newspan
    end

    with_updates_suspended(gl, update = update) do
        gl.nrows += n
        prepend!(gl.rowsizes, rowsizes)
        prepend!(gl.addedrowgaps, addedrowgaps)
    end
end

function prependcols!(gl::GridLayout, n::Int; colsizes=nothing, addedcolgaps=nothing, update = true)

    colsizes = convert_contentsizes(n, colsizes)
    addedcolgaps = convert_gapsizes(n, addedcolgaps, gl.default_colgap)

    foreach(gl.content) do gc
        span = gc.span
        newspan = Span(span.rows, span.cols .+ n)
        gc.span = newspan
    end

    with_updates_suspended(gl, update = update) do
        gl.ncols += n
        prepend!(gl.colsizes, colsizes)
        prepend!(gl.addedcolgaps, addedcolgaps)
    end
end

"""
    insertrows!(gl::GridLayout, at::Int, n::Int; rowsizes=nothing, addedrowgaps=nothing)

Insert `n` rows at row `at` into `GridLayout` `gl`. The new row sizes and row gaps can be
optionally set with `rowsizes` and `addedrowgaps` keywords.
Objects spanning from at least `at-1` up to or beyond `at` are getting extended to span
over the new rows. Objects from `at` and beyond are pushed back, objects before `at` are
unaffected.
"""
function insertrows!(gl::GridLayout, at::Int, n::Int; rowsizes=nothing, addedrowgaps=nothing)

    if !(1 <= at <= nrows(gl))
        error("Invalid row insertion at row $at. GridLayout has $(nrows(gl)) rows.")
    end

    rowsizes = convert_contentsizes(n, rowsizes)
    addedrowgaps = convert_gapsizes(n, addedrowgaps, gl.default_rowgap)

    foreach(gl.content) do gc
        span = gc.span
        rows = span.rows
        newrows = if rows.start < at <= rows.stop
            rows.start : rows.stop + n
        elseif rows.stop < at
            rows
        elseif rows.start >= at
            rows .+ n
        end
        newspan = Span(newrows, span.cols)
        gc.span = newspan
    end

    with_updates_suspended(gl) do
        gl.nrows += n
        splice!(gl.rowsizes, at:at-1, rowsizes)
        splice!(gl.addedrowgaps, at:at-1, addedrowgaps)
    end
end

"""
    insertcols!(gl::GridLayout, at::Int, n::Int; colsizes=nothing, addedcolgaps=nothing)

Insert `n` columns at column `at` into `GridLayout` `gl`. The new column sizes and column gaps can be
optionally set with `colsizes` and `addedcolgaps` keywords.
Objects spanning from at least `at-1` up to or beyond `at` are getting extended to span
over the new columns. Objects from `at` and beyond are pushed back, objects before `at` are
unaffected.
"""
function insertcols!(gl::GridLayout, at::Int, n::Int; colsizes=nothing, addedcolgaps=nothing)

    if !(1 <= at <= ncols(gl))
        error("Invalid column insertion at column $at. GridLayout has $(ncols(gl)) columns.")
    end

    colsizes = convert_contentsizes(n, colsizes)
    addedcolgaps = convert_gapsizes(n, addedcolgaps, gl.default_colgap)

    foreach(gl.content) do gc
        span = gc.span
        cols = span.cols
        newcols = if cols.start < at <= cols.stop
            cols.start : cols.stop + n
        elseif cols.stop < at
            cols
        elseif cols.start >= at
            cols .+ n
        end
        newspan = Span(span.rows, newcols)
        gc.span = newspan
    end

    with_updates_suspended(gl) do
        gl.ncols += n
        splice!(gl.colsizes, at:at-1, colsizes)
        splice!(gl.addedcolgaps, at:at-1, addedcolgaps)
    end
end

function deleterow!(gl::GridLayout, irow::Int)
    if !(1 <= irow <= gl.nrows)
        error("Row $irow does not exist.")
    end

    if gl.nrows == 1
        error("Can't delete the last row")
    end

    # new_content = GridContent[]
    to_remove = GridContent[]
    for c in gl.content
        rows = c.span.rows
        newrows = if irow in rows
            # range is one shorter now
            rows.start : rows.stop - 1
        elseif irow > rows.stop
            # content before deleted row stays the same
            rows
        else
            # content completely after is moved forward 1 step
            rows .- 1
        end
        if isempty(newrows)
            # the row span was just one row and now zero, remove the element
            push!(to_remove, c)
        else
            c.span = Span(newrows, c.span.cols)
        end
    end

    for c in to_remove
        remove_from_gridlayout!(c)
    end
    # gl.content = new_content
    deleteat!(gl.rowsizes, irow)
    deleteat!(gl.addedrowgaps, irow == 1 ? 1 : irow - 1)
    gl.nrows -= 1
    gl.needs_update[] = true
end

function deletecol!(gl::GridLayout, icol::Int)
    if !(1 <= icol <= gl.ncols)
        error("Col $icol does not exist.")
    end

    if gl.ncols == 1
        error("Can't delete the last col")
    end

    to_remove = GridContent[]
    for c in gl.content
        cols = c.span.cols
        newcols = if icol in cols
            # range is one shorter now
            cols.start : cols.stop - 1
        elseif icol > cols.stop
            # content before deleted col stays the same
            cols
        else
            # content completely after is moved forward 1 step
            cols .- 1
        end
        if isempty(newcols)
            # the col span was just one col and now zero, remove the element
            push!(to_remove, c)
        else
            c.span = Span(c.span.rows, newcols)
        end
    end

    for c in to_remove
        remove_from_gridlayout!(c)
    end
    # gl.content = new_content
    deleteat!(gl.colsizes, icol)
    deleteat!(gl.addedcolgaps, icol == 1 ? 1 : icol - 1)
    gl.ncols -= 1
    gl.needs_update[] = true
end

function Base.isempty(gl::GridLayout, dir::GridDir, i::Int)
    !any(gl.content) do c
        span = dir isa Row ? c.span.rows : c.span.cols
        i in span
    end
end

function trim!(gl::GridLayout)
    irow = 1
    while irow <= gl.nrows && gl.nrows > 1
        if isempty(gl, Row(), irow)
            deleterow!(gl, irow)
        else
            irow += 1
        end
    end

    icol = 1
    while icol <= gl.ncols && gl.ncols > 1
        if isempty(gl, Col(), icol)
            deletecol!(gl, icol)
        else
            icol += 1
        end
    end
end

function gridnest!(gl::GridLayout, rows::Indexables, cols::Indexables)

    newrows, newcols = adjust_rows_cols!(gl, rows, cols)

    subgl = GridLayout(
        length(newrows), length(newcols);
        parent = nothing,
        colsizes = gl.colsizes[newcols],
        rowsizes = gl.rowsizes[newrows],
        addedrowgaps = gl.addedrowgaps[newrows.start:(newrows.stop-1)],
        addedcolgaps = gl.addedcolgaps[newcols.start:(newcols.stop-1)],
    )

    # remove the content from the parent that is completely inside the replacement grid
    subgl.block_updates = true
    i = 1
    while i <= length(gl.content)
        gc = gl.content[i]

        if (gc.span.rows.start >= newrows.start && gc.span.rows.stop <= newrows.stop &&
            gc.span.cols.start >= newcols.start && gc.span.cols.stop <= newcols.stop)

            subgl[gc.span.rows .- (newrows.start - 1), gc.span.cols .- (newcols.start - 1), gc.side] = gc.content
            continue
            # don't advance i because there's one piece of content less in the queue
            # and the next item is in the same position as the old removed one
        end
        i += 1
    end
    subgl.block_updates = false

    gl[newrows, newcols] = subgl

    subgl
end

function Base.show(io::IO, ::MIME"text/plain", gl::GridLayout)

    function spaceindent(str, n, downconnection)
        joinstr = if downconnection
            "\n" * (" " ^ 1) * "┃" * (" " ^ (n-2))
        else
            "\n" * (" " ^ n)
        end
        join(split(str, "\n"), joinstr)
    end

    println(io, "GridLayout[$(gl.nrows), $(gl.ncols)] with $(length(gl.content)) children")

    simplespan(span) = span.start == span.stop ? span.start : span

    for (i, c) in enumerate(gl.content)
        rows = c.span.rows
        cols = c.span.cols
        content = c.content

        connector = i == length(gl.content) ? " ┗━ " : " ┣━ "

        if content isa GridLayout
            downconnection = i < length(gl.content)
            str = spaceindent(repr(MIME"text/plain"(), content), 2, downconnection)
            println(io, connector * "[$(simplespan(rows)), $(simplespan(cols))] $str")
        else
            println(io, connector * "[$(simplespan(rows)), $(simplespan(cols))] $(typeof(content))")
        end
    end
end

function Base.show(io::IO, gl::GridLayout)
    print(io, "GridLayout[$(gl.nrows), $(gl.ncols)] ($(length(gl.content)) children)")
end

function colsize!(gl::GridLayout, i::Int, s::ContentSize)
    if !(1 <= i <= gl.ncols)
        error("Can't set size of invalid column $i.")
    end
    gl.colsizes[i] = s
    gl.needs_update[] = true
end

colsize!(gl::GridLayout, i::Int, s::Real) = colsize!(gl, i, Fixed(s))

function rowsize!(gl::GridLayout, i::Int, s::ContentSize)
    if !(1 <= i <= gl.nrows)
        error("Can't set size of invalid row $i.")
    end
    gl.rowsizes[i] = s
    gl.needs_update[] = true
end

rowsize!(gl::GridLayout, i::Int, s::Real) = rowsize!(gl, i, Fixed(s))

function colgap!(gl::GridLayout, i::Int, s::GapSize)
    if !(1 <= i <= (gl.ncols - 1))
        error("Can't set size of invalid column gap $i.")
    end
    gl.addedcolgaps[i] = s
    gl.needs_update[] = true
end

colgap!(gl::GridLayout, i::Int, s::Real) = colgap!(gl, i, Fixed(s))

function colgap!(gl::GridLayout, s::GapSize)
    gl.addedcolgaps .= Ref(s)
    gl.needs_update[] = true
end

function colgap!(gl::GridLayout, r::Real)
    gl.addedcolgaps .= Ref(Fixed(r))
    gl.needs_update[] = true
end

function rowgap!(gl::GridLayout, i::Int, s::GapSize)
    if !(1 <= i <= (gl.nrows - 1))
        error("Can't set size of invalid row gap $i.")
    end
    gl.addedrowgaps[i] = s
    gl.needs_update[] = true
end

rowgap!(gl::GridLayout, i::Int, s::Real) = rowgap!(gl, i, Fixed(s))

function rowgap!(gl::GridLayout, s::GapSize)
    gl.addedrowgaps .= Ref(s)
    gl.needs_update[] = true
end

function rowgap!(gl::GridLayout, r::Real)
    gl.addedrowgaps .= Ref(Fixed(r))
    gl.needs_update[] = true
end

"""
This function solves a grid layout such that the "important lines" fit exactly
into a given bounding box. This means that the protrusions of all objects inside
the grid are not taken into account. This is needed if the grid is itself placed
inside another grid.
"""
function align_to_bbox!(gl::GridLayout, suggestedbbox::Rect2f)

    # compute the actual bbox for the content given that there might be outside
    # padding that needs to be removed
    alignmode = gl.alignmode[]
    bbox = if alignmode isa Outside
        pad = alignmode.padding
        BBox(
            left(suggestedbbox) + pad.left,
            right(suggestedbbox) - pad.right,
            bottom(suggestedbbox) + pad.bottom,
            top(suggestedbbox) - pad.top)
    elseif alignmode isa Inside
        suggestedbbox
    elseif alignmode isa Mixed
        sides = alignmode.sides
        BBox(
            left(suggestedbbox) + (sides.left isa Float32 ? sides.left : 0f0),
            right(suggestedbbox) - (sides.right isa Float32 ? sides.right : 0f0),
            bottom(suggestedbbox) + (sides.bottom isa Float32 ? sides.bottom : 0f0),
            top(suggestedbbox) - (sides.top isa Float32 ? sides.top : 0f0))
    else
        error("Unknown AlignMode of type $(typeof(alignmode))")
    end

    # first determine how big the protrusions on each side of all columns and rows are
    maxgrid = RowCols(gl.ncols, gl.nrows)
    # go through all the layout objects placed in the grid
    for c in gl.content
        idx_rect = side_indices(c)
        mapsides(idx_rect, maxgrid) do side, idx, grid
            grid[idx] = max(grid[idx], protrusion(c, side))
        end
    end

    # for the outside alignmode
    topprot = maxgrid.tops[1]
    bottomprot = maxgrid.bottoms[end]
    leftprot = maxgrid.lefts[1]
    rightprot = maxgrid.rights[end]

    # compute what size the gaps between rows and columns need to be
    colgaps = maxgrid.lefts[2:end] .+ maxgrid.rights[1:end-1]
    rowgaps = maxgrid.tops[2:end] .+ maxgrid.bottoms[1:end-1]

    # determine the biggest gap
    # using the biggest gap size for all gaps will make the layout more even
    if gl.equalprotrusiongaps[2]
        colgaps = ones(gl.ncols - 1) .* (gl.ncols <= 1 ? 0.0 : maximum(colgaps))
    end
    if gl.equalprotrusiongaps[1]
        rowgaps = ones(gl.nrows - 1) .* (gl.nrows <= 1 ? 0.0 : maximum(rowgaps))
    end

    # determine the vertical and horizontal space needed just for the gaps
    # again, the gaps are what the protrusions stick into, so they are not actually "empty"
    # depending on what sticks out of the plots
    sumcolgaps = (gl.ncols <= 1) ? 0.0 : sum(colgaps)
    sumrowgaps = (gl.nrows <= 1) ? 0.0 : sum(rowgaps)

    # compute what space remains for the inner parts of the plots
    remaininghorizontalspace = if alignmode isa Inside
        width(bbox) - sumcolgaps
    elseif alignmode isa Outside
        width(bbox) - sumcolgaps - leftprot - rightprot
    elseif alignmode isa Mixed
        rightal = getside(alignmode, Right())
        leftal = getside(alignmode, Left())
        width(bbox) - sumcolgaps -
            (isnothing(leftal) ? zero(leftprot) : isa(leftal, Protrusion) ? leftal.p : leftprot) -
            (isnothing(rightal) ? zero(rightprot) : isa(rightal, Protrusion) ? rightal.p : rightprot)
    else
        error("Unknown AlignMode of type $(typeof(alignmode))")
    end

    remainingverticalspace = if alignmode isa Inside
        height(bbox) - sumrowgaps
    elseif alignmode isa Outside
        height(bbox) - sumrowgaps - topprot - bottomprot
    elseif alignmode isa Mixed
        topal = getside(alignmode, Top())
        bottomal = getside(alignmode, Bottom())
        height(bbox) - sumrowgaps -
            (isnothing(bottomal) ? zero(bottomprot) : isa(bottomal, Protrusion) ? bottomal.p : bottomprot) -
            (isnothing(topal) ? zero(topprot) : isa(topal, Protrusion) ? topal.p : topprot)
    else
        error("Unknown AlignMode of type $(typeof(alignmode))")
    end

    # compute how much gap to add, in case e.g. labels are too close together
    # this is given as a fraction of the space used for the inner parts of the plots
    # so far, but maybe this should just be an absolute pixel value so it doesn't change
    # when resizing the window
    addedcolgaps = map(gl.addedcolgaps) do cg
        if cg isa Fixed
            return cg.x
        elseif cg isa Relative
            return cg.x * remaininghorizontalspace
        else
            return 0.0 # for float type inference
        end
    end
    addedrowgaps = map(gl.addedrowgaps) do rg
        if rg isa Fixed
            return rg.x
        elseif rg isa Relative
            return rg.x * remainingverticalspace
        else
            return 0.0 # for float type inference
        end
    end

    # compute the actual space available for the rows and columns (plots without protrusions)
    spaceforcolumns = remaininghorizontalspace - ((gl.ncols <= 1) ? 0.0 : sum(addedcolgaps))
    spaceforrows = remainingverticalspace - ((gl.nrows <= 1) ? 0.0 : sum(addedrowgaps))

    colwidths, rowheights = compute_col_row_sizes(spaceforcolumns, spaceforrows, gl)

    # don't allow smaller widths than 1 px even if it breaks the layout (better than weird glitches)
    colwidths = max.(colwidths, ones(length(colwidths)))
    rowheights = max.(rowheights, ones(length(rowheights)))

    # this is the vertical / horizontal space between the inner lines of all plots
    finalcolgaps = colgaps .+ addedcolgaps
    finalrowgaps = rowgaps .+ addedrowgaps

    # compute the resulting width and height of the gridlayout and compute
    # adjustments for the grid's alignment (this will only matter if the grid is
    # bigger or smaller than the bounding box it occupies)

    gridwidth = sum(colwidths) + sum(finalcolgaps) +
        (alignmode isa Outside ? (leftprot + rightprot) : 0.0)
    gridheight = sum(rowheights) + sum(finalrowgaps) +
        (alignmode isa Outside ? (topprot + bottomprot) : 0.0)

    # compute the x values for all left and right column boundaries
    xleftcols = if alignmode isa Inside
        left(bbox) .+ zcumsum(colwidths[1:end-1]) .+
            zcumsum(finalcolgaps)
    elseif alignmode isa Outside
        left(bbox) .+ zcumsum(colwidths[1:end-1]) .+
            zcumsum(finalcolgaps) .+ leftprot
    elseif alignmode isa Mixed
        leftal = getside(alignmode, Left())
        left(bbox) .+ zcumsum(colwidths[1:end-1]) .+
            zcumsum(finalcolgaps) .+ (isnothing(leftal) ? zero(leftprot) : isa(leftal, Protrusion) ? leftal.p : leftprot)
    else
        error("Unknown AlignMode of type $(typeof(alignmode))")
    end
    xrightcols = xleftcols .+ colwidths

    # compute the y values for all top and bottom row boundaries
    ytoprows = if alignmode isa Inside
        top(bbox) .- zcumsum(rowheights[1:end-1]) .-
            zcumsum(finalrowgaps)
    elseif alignmode isa Outside
        top(bbox) .- zcumsum(rowheights[1:end-1]) .-
            zcumsum(finalrowgaps) .- topprot
    elseif alignmode isa Mixed
        topal = getside(alignmode, Top())
        top(bbox) .- zcumsum(rowheights[1:end-1]) .-
            zcumsum(finalrowgaps) .- (isnothing(topal) ? zero(topprot) : isa(topal, Protrusion) ? topal.p : topprot)
    else
        error("Unknown AlignMode of type $(typeof(alignmode))")
    end
    ybottomrows = ytoprows .- rowheights

    # now we can solve the content thats inside the grid because we know where each
    # column and row is placed, how wide it is, etc.
    # note that what we did at the top was determine the protrusions of all grid content,
    # but we know the protrusions before we know how much space each plot actually has
    # because the protrusions should be static (like tick labels etc don't change size with the plot)

    gridboxes = RowCols(
        xleftcols, xrightcols,
        ytoprows, ybottomrows
    )

    for c in gl.content
        idx_rect = side_indices(c)
        bbox_cell = mapsides(idx_rect, gridboxes) do side, idx, gridside
            gridside[idx]
        end

        solving_bbox = bbox_for_solving_from_side(maxgrid, bbox_cell, idx_rect, c.side)

        suggestedbboxobservable(c.content)[] = solving_bbox
    end

    nothing
end

dirlength(gl::GridLayout, c::Col) = gl.ncols
dirlength(gl::GridLayout, r::Row) = gl.nrows

function dirgaps(gl::GridLayout, dir::GridDir)
    starts = zeros(Float32, dirlength(gl, dir))
    stops = zeros(Float32, dirlength(gl, dir))
    for c in gl.content
        span = getspan(c, dir)
        start = span.start
        stop = span.stop
        starts[start] = max(starts[start], protrusion(c, startside(dir)))
        stops[stop] = max(stops[stop], protrusion(c, stopside(dir)))
    end
    starts, stops
end

dirsizes(gl::GridLayout, c::Col) = gl.colsizes
dirsizes(gl::GridLayout, r::Row) = gl.rowsizes

"""
Determine the size of a grid layout along one of its dimensions.
`Row` measures from bottom to top and `Col` from left to right.
The size is dependent on the alignmode of the grid, `Outside` includes
protrusions and paddings.
"""
function determinedirsize(gl::GridLayout, gdir::GridDir)
    sum_dirsizes = 0

    sizes = dirsizes(gl, gdir)

    for idir in 1:dirlength(gl, gdir)
        # width can only be determined for fixed and auto
        sz = sizes[idir]
        dsize = determinedirsize(idir, gl, gdir)

        if isnothing(dsize)
            # early exit if a colsize can not be determined
            return nothing
        end
        sum_dirsizes += dsize
    end

    dirgapsstart, dirgapsstop = dirgaps(gl, gdir)

    forceequalprotrusiongaps = gl.equalprotrusiongaps[gdir isa Row ? 1 : 2]

    dirgapsizes = if forceequalprotrusiongaps
        innergaps = dirgapsstart[2:end] .+ dirgapsstop[1:end-1]
        m = maximum(innergaps)
        innergaps .= m
    else
        innergaps = dirgapsstart[2:end] .+ dirgapsstop[1:end-1]
    end

    inner_gapsizes = dirlength(gl, gdir) > 1 ? sum(dirgapsizes) : 0

    addeddirgapsizes = gdir isa Row ? gl.addedrowgaps : gl.addedcolgaps

    addeddirgaps = dirlength(gl, gdir) == 1 ? 0 : sum(addeddirgapsizes) do c
        if c isa Fixed
            c.x
        elseif c isa Relative
            error("Auto grid size not implemented with relative gaps")
        end
    end

    inner_size_combined = sum_dirsizes + inner_gapsizes + addeddirgaps
    return if gl.alignmode[] isa Inside
        inner_size_combined
    elseif gl.alignmode[] isa Outside
        paddings = if gdir isa Row
            gl.alignmode[].padding.top + gl.alignmode[].padding.bottom
        else
            gl.alignmode[].padding.left + gl.alignmode[].padding.right
        end
        inner_size_combined + dirgapsstart[1] + dirgapsstop[end] + paddings
    else
        error("Unknown AlignMode of type $(typeof(gl.alignmode[]))")
    end
end

"""
Determine the size of one row or column of a grid layout.
"""
function determinedirsize(idir, gl, dir::GridDir)

    sz = dirsizes(gl, dir)[idir]

    if sz isa Fixed
        # fixed dir size can simply be returned
        return sz.x
    elseif sz isa Relative
        # relative dir size can't be inferred
        return nothing
    elseif sz isa Auto
        # auto dir size can either be determined or not, depending on the
        # trydetermine flag
        !sz.trydetermine && return nothing

        dirsize = nothing
        for c in gl.content
            # content has to be single span to be determinable in size
            singlespanned = getspan(c, dir).start == getspan(c, dir).stop == idir

            # content has to be placed with Inner side, otherwise it's protrusion
            # content
            is_inner = c.side isa Inner

            if singlespanned && is_inner
                s = determinedirsize(c.content, dir, c.side)
                if !isnothing(s)
                    dirsize = isnothing(dirsize) ? s : max(dirsize, s)
                end
            end
        end
        return dirsize
    end
    nothing
end

# a function that iterates over those sizes that belong to a type T
# while enumerating all indices, so that i can be used to index colwidths / rowheights
# and determinedcols / determinedrows
filterenum(f, T::Type, iter) = foreach(f, ((i, value) for (i, value) in enumerate(iter) if value isa T))

function compute_col_row_sizes(spaceforcolumns, spaceforrows, gl)
    # the space for columns and for rows is divided depending on the sizes
    # stored in the grid layout

    # algorithm:

    # 1. get fixed sizes
    # 2. compute relative sizes
    # 3. determine determinable auto sizes
    # 4. compute those aspect sizes that are relative to one of the three above categories
    # 5. at least one side now has to have only undeterminable auto sizes left
    # 6. compute remaining auto sizes for one side
    # 7. compute remaining aspect sizes on other side
    # 8. compute remaining auto sizes on the same side

    colwidths = zeros(gl.ncols)
    rowheights = zeros(gl.nrows)

    determinedcols = zeros(Bool, gl.ncols)
    determinedrows = zeros(Bool, gl.nrows)

    # first fixed sizes
    filterenum(Fixed, gl.colsizes) do (i, fixed)
        colwidths[i] = fixed.x
        determinedcols[i] = true
    end
    filterenum(Fixed, gl.rowsizes) do (i, fixed)
        rowheights[i] = fixed.x
        determinedrows[i] = true
    end

    # then relative sizes
    filterenum(Relative, gl.colsizes) do (i, relative)
        colwidths[i] = relative.x * spaceforcolumns
        determinedcols[i] = true
    end
    filterenum(Relative, gl.rowsizes) do (i, relative)
        rowheights[i] = relative.x * spaceforrows
        determinedrows[i] = true
    end

    # then determinable auto sizes
    filterenum(Auto, gl.colsizes) do (i, auto)
        size = determinedirsize(i, gl, Col())
        if !isnothing(size)
            colwidths[i] = size
            determinedcols[i] = true
        end
    end
    filterenum(Auto, gl.rowsizes) do (i, auto)
        size = determinedirsize(i, gl, Row())
        if !isnothing(size)
            rowheights[i] = size
            determinedrows[i] = true
        end
    end

    # now aspect sizes that refer to already determined counterparts
    filterenum(Aspect, gl.colsizes) do (i, aspect)
        if determinedrows[aspect.index]
            colwidths[i] = aspect.ratio * rowheights[aspect.index]
            determinedcols[i] = true
        end
    end
    filterenum(Aspect, gl.rowsizes) do (i, aspect)
        if determinedcols[aspect.index]
            rowheights[i] = aspect.ratio * colwidths[aspect.index]
            determinedrows[i] = true
        end
    end

    remaining_colspace = spaceforcolumns - sum(colwidths)
    remaining_rowspace = spaceforrows - sum(rowheights)

    # if we have aspect sizes left on one side, they can only be determined
    # if the other side has only undeterminable autos left
    n_col_aspects_left = sum(enumerate(gl.colsizes)) do (i, size)
        (size isa Aspect) && (determinedcols[i] == false)
    end
    n_row_aspects_left = sum(enumerate(gl.rowsizes)) do (i, size)
        (size isa Aspect) && (determinedrows[i] == false)
    end

    n_col_autos_left = sum(enumerate(gl.colsizes)) do (i, size)
        (size isa Auto) && (determinedcols[i] == false)
    end
    n_row_autos_left = sum(enumerate(gl.rowsizes)) do (i, size)
        (size isa Auto) && (determinedrows[i] == false)
    end

    if n_col_aspects_left == 0
        let
            indices = Int[]
            ratios = Float64[]
            i_ratios = filterenum(Auto, gl.colsizes) do (i, auto)
                if determinedcols[i] == false
                    push!(indices, i)
                    push!(ratios, auto.ratio)
                end
            end
            sumratios = sum(ratios)
            for (i, ratio) in zip(indices, ratios)
                colwidths[i] = ratio / sumratios * remaining_colspace
                determinedcols[i] = true
            end
        end
    end

    if n_row_aspects_left == 0
        let
            indices = Int[]
            ratios = Float64[]
            i_ratios = filterenum(Auto, gl.rowsizes) do (i, auto)
                if determinedrows[i] == false
                    push!(indices, i)
                    push!(ratios, auto.ratio)
                end
            end
            sumratios = sum(ratios)
            for (i, ratio) in zip(indices, ratios)
                rowheights[i] = ratio / sumratios * remaining_rowspace
                determinedrows[i] = true
            end
        end
    end

    # now if either columns or rows had no aspects left, they should have all sizes determined
    # we run over the aspects again
    filterenum(Aspect, gl.colsizes) do (i, aspect)
        if determinedrows[aspect.index]
            colwidths[i] = aspect.ratio * rowheights[aspect.index]
            determinedcols[i] = true
        else
            error("Column $i was given an Aspect size relative to row $(aspect.index). This row's size could not be determined in time, therefore the layouting algorithm failed. This probably happened because you used an Aspect row and column size at the same time, which couldn't both be resolved.")
        end
    end
    filterenum(Aspect, gl.rowsizes) do (i, aspect)
        if determinedcols[aspect.index]
            rowheights[i] = aspect.ratio * colwidths[aspect.index]
            determinedrows[i] = true
        else
            error("Row $i was given an Aspect size relative to column $(aspect.index). This column's size could not be determined in time, therefore the layouting algorithm failed. This probably happened because you used an Aspect row and column size at the same time, which couldn't both be resolved.")
        end
    end

    # if we haven't errored yet, all aspect sizes are done
    # one more pass over the undetermined autos is all that's needed

    remaining_colspace = spaceforcolumns - sum(colwidths)
    remaining_rowspace = spaceforrows - sum(rowheights)

    let
        indices = Int[]
        ratios = Float64[]
        i_ratios = filterenum(Auto, gl.colsizes) do (i, auto)
            if determinedcols[i] == false
                push!(indices, i)
                push!(ratios, auto.ratio)
            end
        end
        sumratios = sum(ratios)
        for (i, ratio) in zip(indices, ratios)
            colwidths[i] = ratio / sumratios * remaining_colspace
            determinedcols[i] = true
        end
    end

    let
        indices = Int[]
        ratios = Float64[]
        i_ratios = filterenum(Auto, gl.rowsizes) do (i, auto)
            if determinedrows[i] == false
                push!(indices, i)
                push!(ratios, auto.ratio)
            end
        end
        sumratios = sum(ratios)
        for (i, ratio) in zip(indices, ratios)
            rowheights[i] = ratio / sumratios * remaining_rowspace
            determinedrows[i] = true
        end
    end

    # now all columns and rows should have their sizes
    ncols_undetermined = sum(.!determinedcols)
    nrows_undetermined = sum(.!determinedrows)

    if ncols_undetermined > 0
        error("After a non-erroring layouting pass, the number of undetermined columns is $ncols_undetermined. This must be a bug.")
    end
    if nrows_undetermined > 0
        error("After a non-erroring layouting pass, the number of undetermined rows is $nrows_undetermined. This must be a bug.")
    end

    colwidths, rowheights
end

function Base.setindex!(g::GridLayout, content, rows::Indexables, cols::Indexables, side::Side = Inner())
    add_content!(g, content, rows, cols, side)
    content
end

function Base.setindex!(g::GridLayout, content_array::AbstractArray{T, 2}) where T
    rowrange = 1:size(content_array, 1)
    colrange = 1:size(content_array, 2)
    g[rowrange, colrange] = content_array
end

function Base.setindex!(g::GridLayout, content_array::AbstractArray{T, 1}) where T
    error("""
        You can only assign a one-dimensional content AbstractArray if you also specify the direction in the layout.
        Valid options are :h for horizontal and :v for vertical.
        Example:
            layout[:h] = contentvector
    """)
end

function Base.setindex!(g::GridLayout, content_array::AbstractArray{T, 1}, h_or_v::Symbol) where T
    if h_or_v == :h
        g[1, 1:length(content_array)] = content_array
    elseif h_or_v == :v
        g[1:length(content_array), 1] = content_array
    else
        error("""
            Invalid direction specifier $h_or_v.
            Valid options are :h for horizontal and :v for vertical.
        """)
    end
end

function Base.setindex!(g::GridLayout, content_array::AbstractArray, rows::Indexables, cols::Indexables)

    rows, cols = to_ranges(g, rows, cols)

    if rows.start < 1
        error("Can't prepend rows using array syntax so far, start row $(rows.start) is smaller than 1.")
    end
    if cols.start < 1
        error("Can't prepend columns using array syntax so far, start column $(cols.start) is smaller than 1.")
    end

    nrows = length(rows)
    ncols = length(cols)
    ncells = nrows * ncols

    if ndims(content_array) == 2
        if size(content_array) != (nrows, ncols)
            error("Content array size is size $(size(content_array)) for $nrows rows and $ncols cols")
        end
        # put the array content into the grid layout in order
        for (i, r) in enumerate(rows), (j, c) in enumerate(cols)
            g[r, c] = content_array[i, j]
        end
    elseif ndims(content_array) == 1
        if length(content_array) != nrows * ncols
            error("Content array size is length $(length(content_array)) for $nrows * $ncols cells")
        end
        # put the content in the layout along columns first, because that is more
        # intuitive
        for (i, (c, r)) in enumerate(Iterators.product(cols, rows))
            g[r, c] = content_array[i]
        end
    else
        error("Can't assign a content array with $(ndims(content_array)) dimensions, only 1 or 2.")
    end
    content_array
end

function GridContent(content::T, span::Span, side::Side) where T
    needs_update = Observable(false)
    # connect the correct observables
    protrusions_handle = on(protrusionsobservable(content)) do p
        needs_update[] = true
    end
    reportedsize_handle = on(reportedsizeobservable(content)) do c
        needs_update[] = true
    end
    GridContent{GridLayout, T}(nothing, content, span, side, needs_update,
        protrusions_handle, reportedsize_handle)
end

function add_content!(g::GridLayout, content, rows, cols, side::Side)
    # update = false because update is called in add_to_gridlayout! anyway
    rows, cols = adjust_rows_cols!(g, rows, cols; update = false)

    gc = if !isnothing(gridcontent(content))
        # take the existing gridcontent, remove it from its gridlayout if it has one,
        # and modify it with the new span and side
        gridc = gridcontent(content)
        remove_from_gridlayout!(gridc)
        gridc.span = Span(rows, cols)
        gridc.side = side
        gridc
    else
        # make a new one if none existed
        GridContent(content, Span(rows, cols), side)
    end

    layoutobservables(content).gridcontent = gc

    connect_layoutobservables!(gc)

    add_to_gridlayout!(g, gc)
end

function Base.lastindex(g::GridLayout, d)
    if d == 1
        g.nrows
    elseif d == 2
        g.ncols
    else
        error("A grid only has two dimensions, you're indexing dimension $d.")
    end
end

function GridPosition(g::GridLayout, rows::Indexables, cols::Indexables, side = Inner())
    span = Span(to_ranges(g, rows, cols)...)
    GridPosition(g, span, side)
end

function Base.getindex(g::GridLayout, rows::Indexables, cols::Indexables, side = Inner())
    GridPosition(g, rows, cols, side)
end

function Base.setindex!(gp::GridPosition, element)
    gp.layout[gp.span.rows, gp.span.cols, gp.side] = element
end

function Base.setindex!(gp::GridPosition, element, rows, cols, side = Inner())
    layout = get_layout_at!(gp, createmissing = true)
    layout[rows, cols, side] = element
    element
end

ncols(g::GridLayout) = g.ncols
nrows(g::GridLayout) = g.nrows
Base.size(g::GridLayout) = (nrows(g), ncols(g))

Base.in(span1::Span, span2::Span) = span1.rows.start >= span2.rows.start &&
    span1.rows.stop <= span2.rows.stop &&
    span1.cols.start >= span2.cols.start &&
    span1.cols.stop <= span2.cols.stop

"""
    contents(gp::GridPosition; exact::Bool = false)

Retrieve all objects placed in the `GridLayout` at the `Span` and `Side` stored
in the `GridPosition` `gp`. If `exact == true`, elements are only included
if they match the `Span` exactly, otherwise they can also be contained within the spanned layout area.
"""
function contents(gp::GridPosition; exact::Bool = false)
    contents = []
    for c in gp.layout.content
        if exact
            if c.span == gp.span && c.side == gp.side
                push!(contents, c.content)
            end
        else
            if c.span in gp.span && c.side == gp.side
                push!(contents, c.content)
            end
        end
    end
    contents
end

"""
    contents(g::GridLayout)

Retrieve all objects placed in the `GridLayout` `g`, in the order they are stored, extracted from
their containing `GridContent`s.
"""
function contents(g::GridLayout)
    map(g.content) do gc
        gc.content
    end
end

function Base.getindex(gp::Union{GridPosition, GridSubposition}, rows, cols, side = Inner())
    GridSubposition(gp, rows, cols, side)
end

function Base.setindex!(parent::GridSubposition, obj,
    rows, cols, side = GridLayoutBase.Inner())
    layout = get_layout_at!(parent, createmissing = true)
    layout[rows, cols, side] = obj
    obj
end

function Base.setindex!(parent::GridSubposition, obj)
    layout = get_layout_at!(parent.parent, createmissing = true)
    layout[parent.rows, parent.cols, parent.side] = obj
    obj
end

function get_layout_at!(gp::GridPosition; createmissing = false)
    c = contents(gp, exact = true)
    layouts = filter(x -> x isa GridLayoutBase.GridLayout, c)
    if isempty(layouts)
        if createmissing
            return gp[] = GridLayoutBase.GridLayout()
        else
            error("No layout found but `createmissing` is false.")
        end
    elseif length(layouts) == 1
        return first(layouts)
    else
        error("Found more than zero or one GridLayouts at $gp")
    end
end

function get_layout_at!(gsp::GridSubposition; createmissing = false)
    layout = get_layout_at!(gsp.parent; createmissing = createmissing)
    gp = layout[gsp.rows, gsp.cols, gsp.side]
    get_layout_at!(gp, createmissing = createmissing)
end

function contents(g::GridSubposition; exact = false)
    layout = get_layout_at!(g.parent, createmissing = false)
    contents(layout[g.rows, g.cols, g.side], exact = exact)
end

function content(g::Union{GridPosition,GridSubposition})
    cs = contents(g, exact = true)
    if length(cs) == 1
        return cs[1]
    else
        error("There is not exactly one object at the given GridPosition")
    end
end

function parent(g::GridLayout)
    g.parent
end

function top_parent(g::GridLayout)
    top_parent(parent(g))
end

top_parent(x) = x

function top_parent_grid(g::GridLayout)
    p = parent(g)
    if p isa GridLayout
        top_parent_grid(p)
    else
        g
    end
end


text:
#-------------------------------------------------------------------------# StatLearn
type StatLearn{
        A <: Algorithm,
        M <: Model,
        P <: Penalty,
        W <: StochasticWeight
    } <: OnlineStat{XYInput}

    β0::Float64     # intercepts
    β::VecF         # coefficients
    intercept::Bool # should β0 be estimated?
    η::Float64      # constant part of learning rate

    # Storage
    H0::Float64
    G0::Float64
    H::VecF
    G::VecF

    algorithm::A    # determines how updates work
    model::M        # model definition
    penalty::P      # type of penalty
    weight::W       # Weight, may not get used, depending on algorithm
end
function _StatLearn(p::Integer, wgt::Weight;
        model::Model = L2Regression(),
        eta::Real = 1.0,
        penalty::Penalty = NoPenalty(),
        algorithm::Algorithm = SGD(),
        intercept::Bool = true,
    )
    o = StatLearn(0.0, zeros(p), intercept, Float64(eta), _ϵ, _ϵ,
        _ϵ * ones(p), _ϵ * ones(p), algorithm, model, penalty, wgt)
    o
end
function StatLearn(p::Integer, args...; kw...)
    wgt = LearningRate()
    mod = LinearRegression()
    alg = SGD()
    pen = NoPenalty()
    for arg in args
        T = typeof(arg)
        if T <: Weight
            wgt = arg
        elseif T <: Model
            mod = arg
        elseif T <: Algorithm
            alg = arg
        elseif T <: Penalty
            pen = arg
        end
    end
    _StatLearn(p, wgt; model = mod, algorithm = alg, penalty = pen, kw...)
end
function StatLearn(x::AMat, y::AVec, args...; kw...)
    o = StatLearn(size(x, 2), args...;kw...)
    fit!(o, x, y)
end
function StatLearn(x::AMat, y::AVec, b::Integer, args...; kw...)
    o = StatLearn(size(x, 2), args...; kw...)
    fit!(o, x, y, b)
    o
end

function Base.show(io::IO, o::StatLearn)
    printheader(io, "StatLearn")
    o.intercept && print_item(io, "Intercept", o.β0)
    print_item(io, "β", o.β)
    print_item(io, "Model", o.model)
    print_item(io, "Penalty", o.penalty)
    print_item(io, "Algorithm", o.algorithm)
    print_item(io, "η", o.η)
    print_item(io, "Weight", typeof(o.weight))
    print_item(io, "Nobs", nobs(o))
end

coef(o::StatLearn) = o.intercept ? vcat(o.β0, o.β) : o.β
predict(o::StatLearn, x) = predict(o.model, xβ(o, x))
classify(o::StatLearn, x) = classify(o.model, xβ(o, x))

xβ(o::StatLearn, x::AVec) = o.β0 + dot(o.β, x)
xβ(o::StatLearn, x::AMat) = o.β0 + x * o.β

loss(o::StatLearn, x::AVec, y::Real) = loss(o.model, y, xβ(o, x))
loss(o::StatLearn, x::AMat, y::AVec) = loss(o.model, y, xβ(o, x))
cost(o::StatLearn, x::AVec, y::Real) = loss(o.model, y, xβ(o, x)) + penalty(o.penalty, o.β)
cost(o::StatLearn, x::AMat, y::AVec) = loss(o.model, y, xβ(o, x)) + penalty(o.penalty, o.β)

#-------------------------------------------------------------------------------# SGD
immutable SGD <: Algorithm end
updateβ0!(o::StatLearn{SGD}, γ, ηγ, g, ηγg) = (o.β0 -= ηγg)
function updateβ!(o::StatLearn{SGD}, P, β, j, γ, ηγ, gx, ηγgx)
    @inbounds β[j] -= ηγ * (gx + deriv(P, β[j]))
end

#---------------------------------------------------------------------------# AdaGrad
immutable AdaGrad <: Algorithm end
function updateβ0!(o::StatLearn{AdaGrad}, γ, ηγ, g, ηγg)
    o.H0 = smooth(o.H0, g * g, 1 / nups(o.weight))
    o.β0 -= ηγg / sqrt(o.H0)
end
function updateβ!(o::StatLearn{AdaGrad}, P, β, j, γ, ηγ, gx, ηγgx)
    @inbounds o.H[j] = smooth(o.H[j], gx * gx, 1 / nups(o.weight))
    @inbounds step = ηγ / sqrt(o.H[j])
    @inbounds β[j] = prox(P, β[j] - step * gx, step)
end

#--------------------------------------------------------------------------# AdaGrad2
immutable AdaGrad2 <: Algorithm end
function updateβ0!(o::StatLearn{AdaGrad2}, γ, ηγ, g, ηγg)
    o.H0 = smooth(o.H0, g * g, γ)
    o.β0 -= ηγg / sqrt(o.H0)
end
function updateβ!(o::StatLearn{AdaGrad2}, P, β, j, γ, ηγ, gx, ηγgx)
    @inbounds o.H[j] = smooth(o.H[j], gx * gx, γ)
    @inbounds step = ηγ / sqrt(o.H[j])
    @inbounds β[j] = prox(P, β[j] - step * gx, step)
end

#--------------------------------------------------------------------------# AdaDelta
immutable AdaDelta <: Algorithm
    ρ::Float64
    AdaDelta(ρ::Real = .001) = new(ρ)
end
function updateβ0!(o::StatLearn{AdaDelta}, γ, ηγ, g, ηγg)
    o.H0 = smooth(o.H0, g * g, o.algorithm.ρ)
    Δ = sqrt(o.G0 / o.H0) * g
    o.β0 -= Δ
    o.G0 = smooth(o.G0, Δ * Δ, o.algorithm.ρ)
end
function updateβ!(o::StatLearn{AdaDelta}, P, β, j, γ, ηγ, gx, ηγgx)
    @inbounds o.H[j] = smooth(o.H[j], gx * gx, o.algorithm.ρ)
    @inbounds step = sqrt(o.G[j] / o.H[j])
    @inbounds β[j] = prox(P, β[j] - step * gx, step)
    Δ = step * gx
    @inbounds o.G[j] = smooth(o.G[j], Δ * Δ, o.algorithm.ρ)
end

#------------------------------------------------------------------------------# ADAM
immutable ADAM <: Algorithm
    m1::Float64
    m2::Float64
    ADAM(m1::Real = .01, m2::Real = .01) = new(m1, m2)
end
function updateβ0!(o::StatLearn{ADAM}, γ, ηγ, g, ηγg)
    m1, m2, nups = o.algorithm.m1, o.algorithm.m2, o.weight.nups
    ratio = sqrt(1.0 - m2 ^ nups) / (1.0 - m1 ^ nups)
    o.H0 = smooth(o.H0, g, m1)
    o.G0 = smooth(o.G0, g * g, m2)
    o.β0 -= ηγ * ratio * o.H0 / sqrt(o.G0)
end
function updateβ!(o::StatLearn{ADAM}, P, β, j, γ, ηγ, gx, ηγgx)
    m1, m2, nups = o.algorithm.m1, o.algorithm.m2, o.weight.nups
    ratio = sqrt(1.0 - m2 ^ nups) / (1.0 - m1 ^ nups)
    @inbounds o.H[j] = smooth(o.H[j], gx, m1)
    @inbounds o.G[j] = smooth(o.G[j], gx * gx, m2)
    @inbounds step = ηγ * ratio / sqrt(o.G[j])
    @inbounds β[j] = prox(P, β[j] - step * o.H[j], step)
end

#--------------------------------------------------------------------------# Momentum
immutable Momentum <: Algorithm
    α::Float64
    Momentum(α::Real = .1) = new(α)
end
function updateβ0!(o::StatLearn{Momentum}, γ, ηγ, g, ηγg)
    o.H0 = smooth(o.H0, g, o.algorithm.α)
    o.β0 -= ηγ * o.H0
end
function updateβ!(o::StatLearn{Momentum}, P, β, j, γ, ηγ, gx, ηγgx)
    @inbounds o.H[j] = smooth(o.H[j], gx, o.algorithm.α)
    @inbounds β[j] -= ηγ * o.H[j]
    @inbounds β[j] -= ηγ * (gx + deriv(P, β[j]))
end

# #-------------------------------------------------------------------------------# OMM
# # Apparently this is equivalent to SGD...
# immutable OMM <: Algorithm end
# function updateβ0!(o::StatLearn{OMM}, γ, ηγ, g, ηγg)
#     o.H0 = smooth(o.H0, g, γ)
#     o.G0 = smooth(o.G0, o.β0, γ)
#     o.β0 = o.G0 - o.η * o.H0
# end
# function updateβ!(o::StatLearn{OMM}, β, j, γ, ηγ, gx, ηγgx)
#     @inbounds o.H[j] = smooth(o.H[j], gx, γ)
#     @inbounds o.G[j] = smooth(o.G[j], o.β[j], γ)
#     @inbounds β[j] = o.G[j] - o.η * o.H[j]
# end

#---------------------------------------------------------------------------# fitting
function _fit!{T <: Real}(o::StatLearn, x::AVec{T}, y::Real, γ::Float64)
    η, β, P = o.η, o.β, o.penalty
    ηγ = η * γ
    xb = dot(x, β) + o.β0
    g = lossderiv(o.model, y, xb)
    ηγg = ηγ * g
    if o.intercept
        updateβ0!(o, γ, ηγ, g, ηγg)
    end
    for j in eachindex(β)
        gx = g * x[j]
        ηγgx = ηγ * gx
        updateβ!(o, P, β, j, γ, ηγ, gx, ηγgx)
    end
    o
end

function _fitbatch!{T<:Real, S<:Real}(o::StatLearn, x::AMat{T}, y::AVec{S}, γ::Float64)
    η, β, P = o.η, o.β, o.penalty
    ηγ = η * γ
    xb = x * β
    gvec = zeros(size(x, 1))
    for i in eachindex(gvec)
        @inbounds gvec[i] = lossderiv(o.model, y[i], xb[i])
    end
    if o.intercept
        g = mean(gvec)
        ηγg = ηγ * g
        updateβ0!(o, γ, ηγ, g, ηγg)
    end
    for j in eachindex(β)
        gx = batch_gx(view(x, :, j), gvec)
        ηγgx = ηγ * gx
        updateβ!(o, P, β, j, γ, ηγ, gx, ηγgx)
    end
    o
end
function batch_gx(xj::AVec, g::AVec)
    v = 0.0
    n = length(xj)
    for i in eachindex(xj)
        @inbounds v += xj[i] * g[i]
    end
    v / n
end

"""
`fit!` and update the tuning parameter for a `StatLearn` object.

`cvfit!(o::StatLearn, x, y, xtest, ytest)`

For regularization parameter `λ` and next weight `γ`, the updated `λ` will be one of:

- `λ + γ`
- `λ`
- `max(λ - γ, 0)`

The choice depends on which results in the smallest loss on the test set `(xtest, ytest)`.
"""
function cvfit!(o::StatLearn, x, y, xtest, ytest, arg = 1)
    @assert typeof(o.penalty) != NoPenalty
    γ = weight(o, size(x, 1))
    h = copy(o); h.penalty = typeof(h.penalty)(o.penalty.λ + γ)
    l = copy(o); l.penalty = typeof(l.penalty)(max(o.penalty.λ - γ, 0.0))
    fit!(h, x, y, arg)
    fit!(o, x, y, arg)
    fit!(l, x, y, arg)
    loss_h = loss(h, xtest, ytest)
    loss_o = loss(o, xtest, ytest)
    loss_l = loss(l, xtest, ytest)
    lossmin = minimum([loss_h, loss_o, loss_l])
    loss_l == lossmin && _copy!(o, l)
    loss_h == lossmin && _copy!(o, h)
    o
end
function _copy!(o::StatLearn, o2::StatLearn)
    o.β0 = o2.β0
    o.β[:] = o2.β
    o.penalty = o2.penalty
end


text:
@with_kw mutable struct MergingNeighborsOverlay <: SceneOverlay
    target_id::Int64 = 0
    env::MergingEnvironment = MergingEnvironment()
    color_fore::Colorant = colorant"blue"
    color_main::Colorant = colorant"green"
    line_width::Float64 = 0.5
    textparams::TextParams = TextParams(y_start=20, size=20)
end

function AutoViz.render!(rendermodel::RenderModel, overlay::MergingNeighborsOverlay,
                 scene::Scene, roadway::Roadway)
    textparams = overlay.textparams
    yₒ = textparams.y_start
    Δy = textparams.y_jump
    vehicle_index = findfirst(overlay.target_id, scene)
    veh_ego = scene[vehicle_index]
    fore, merge, fore_main, rear_main = get_neighbors(overlay.env, scene, overlay.target_id)
    if fore.ind != nothing
        veh_oth = scene[fore.ind]
        A = get_front(veh_ego)
        B = get_rear(veh_oth)
        add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_fore, overlay.line_width))
        # drawtext(@sprintf("d fore:   %10.3f", fore.Δs), yₒ + 0*Δy, rendermodel, textparams)
    end
    if fore_main.ind != nothing
        veh_oth = scene[fore_main.ind]
        A = get_front(veh_ego)
        B = get_rear(veh_oth)
        add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_main, overlay.line_width))
        # drawtext(@sprintf("d fore main: %10.3f", fore_main.Δs), yₒ + 1*Δy, rendermodel, textparams)
    end
    if rear_main.ind != nothing
        veh_oth = scene[rear_main.ind]
        A = get_front(veh_ego)
        B = get_rear(veh_oth)
        add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_main, overlay.line_width))
        # drawtext(@sprintf("d rear main:  %10.3f", rear_main.Δs), yₒ + 2*Δy, rendermodel, textparams)
    end
    if merge.ind != nothing
        veh_oth = scene[merge.ind]
        A = get_front(veh_ego)
        B = get_rear(veh_oth)
        add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_main, overlay.line_width))
        # drawtext(@sprintf("d merge:  %10.3f", merge.Δs), yₒ + 3*Δy, rendermodel, textparams)
    end
end

@with_kw mutable struct DistToMergeOverlay <: SceneOverlay
    target_id::Int64 = 0
    env::MergingEnvironment = MergingEnvironment()
    color_point::Colorant = colorant"red"
    color_line::Colorant = colorant"black"
    point_radius::Float64 = 0.5
    line_width::Float64 = 0.5
    textparams::TextParams = TextParams(y_start=50)
end

function AutoViz.render!(rendermodel::RenderModel, overlay::DistToMergeOverlay,
                 scene::Scene, roadway::Roadway)
    # display merge point
    mp = overlay.env.merge_point
    add_instruction!(rendermodel, render_circle, (mp.x, mp.y, overlay.point_radius, overlay.color_point))

    # display line to merge point + info
    vehind = findfirst(overlay.target_id, scene)
    if vehind != nothing
        veh = scene[vehind]
        # plot line
        A = get_front(veh)
        # add_instruction!(rendermodel, render_dashed_line, ([A.x mp.x; A.y mp.y], overlay.color_line, overlay.line_width))
        dm = dist_to_merge(overlay.env, veh)
        ttm = time_to_merge(overlay.env, veh)
        disp_point = A + VecSE2(-10.0, 10.0, 0.0)
        if overlay.target_id == 1
            disp_point = A + VecSE2(-10.0, -10.0, 0.0)
        end
        text_overlay = TextOverlay(text=[@sprintf("dist=%3.2f m",dm),
                                         @sprintf("vel=%3.2f m/s", veh.state.v),
                                         @sprintf("ttm=%2.2f s", ttm)],
                                   pos=disp_point,
                                   line_spacing = 0.15,
                                   font_size = 20,
                                   incameraframe=true)
        # render!(rendermodel, text_overlay, scene, roadway)
    end
end

@with_kw mutable struct MaskingOverlay <: SceneOverlay
    mdp::GenerativeMergingMDP = GenerativeMergingMDP()
    acc::Float64 = 0.0 # current acceleration
    textparams::TextParams = TextParams(size=20, x = 300, y_start=110)
end

function AutoViz.render!(rendermodel::RenderModel, overlay::MaskingOverlay,
                 scene::Scene, roadway::Roadway)
        ss = AugScene(scene, (acc=overlay.acc,))
        acts = safe_actions(overlay.mdp, ss)
        textparams = overlay.textparams
        drawtext(@sprintf("Available Actions %s", acts), textparams.y_start, rendermodel, textparams)
        # drawtext(@sprintf("Actions %s", acts), textparams.y_start, rendermodel, textparams)

end

@with_kw mutable struct CooperativeIDMOverlay <: SceneOverlay
    model::CooperativeIDM = CooperativeIDM()
    targetid::Int64 = 0
    textparams::TextParams = TextParams(size = 20, x = 700, y_start=30)
end

function AutoViz.render!(rendermodel::RenderModel, overlay::CooperativeIDMOverlay,
    scene::Scene, roadway::Roadway)
    observe!(overlay.model, scene, roadway, overlay.targetid)
    veh = get_by_id(scene, overlay.targetid)
    mergeveh = get_by_id(scene, EGO_ID)
    textparams = overlay.textparams
    yₒ = textparams.y_start
    Δy = textparams.y_jump
    drawtext(@sprintf("Driver Model Info %1d:", overlay.targetid), yₒ, rendermodel, textparams)
    # drawtext(@sprintf("TTM threshold: %1.1f", overlay.model.ttm_threshold), yₒ + Δy, rendermodel, textparams)
    drawtext(@sprintf("cooperation level: %1.1f", overlay.model.c), yₒ + Δy, rendermodel, textparams)
    drawtext(@sprintf("act merge: %1.1f m/s^2", overlay.model.a_merge), yₒ + 2*Δy, rendermodel, textparams)
    drawtext(@sprintf("act idm:   %1.1f m/s^2", overlay.model.a_idm), yₒ + 3*Δy, rendermodel, textparams)
    drawtext(@sprintf("act final: %1.1f m/s^2", overlay.model.a), yₒ + 4*Δy, rendermodel, textparams)
    drawtext(@sprintf("dist at merge: %2.2f", overlay.model.dist_at_merge), yₒ + 5*Δy, rendermodel, textparams)
    drawtext(@sprintf("desired gap: %2.2f", overlay.model.s_des), yₒ + 6*Δy, rendermodel, textparams)
    ego_ttm = overlay.model.ego_ttm
    veh_ttm = overlay.model.veh_ttm
    # @show !(ego_ttm < 0.0 || ego_ttm = Inf || ego_ttm <= veh_ttm || veh_ttm == Inf)
    # drawtext(@sprintf("considering merge: %s", !(ego_ttm < 0.0 || ego_ttm == Inf || ego_ttm < veh_ttm || veh_ttm == Inf || mergeveh.state.posF.s < overlay.model.fov)), yₒ + 7*Δy, rendermodel, textparams)
    drawtext(@sprintf("considering merge: %s", overlay.model.consider_merge), yₒ + 7*Δy, rendermodel, textparams)
    drawtext(@sprintf("driver ttm: %2.1f", overlay.model.ego_ttm), yₒ + 8*Δy, rendermodel, textparams)
    drawtext(@sprintf("merge ttm: %2.1f", overlay.model.veh_ttm), yₒ + 9*Δy, rendermodel, textparams)
    drawtext(@sprintf("front car %s", overlay.model.front_car), yₒ + 10*Δy, rendermodel, textparams)
end

function get_car_type_colors(scene::Scene, models::Dict{Int64, DriverModel};
                             cooperative::Colorant = RGBA(0.0, 1., 0., 1.),
                             aggressive::Colorant = RGBA(1.0,0.0,0.498, 1.),
                             cooperative_slow::Colorant = RGBA(0.8, 0.9, 0., 1.),
                             aggressive_slow::Colorant = RGBA(0.7,0.0,0.9, 1) )
    color_dict = Dict{Int64, Colorant}()
    for veh in scene
        if veh.id == 1
            color_dict[1] = AutoViz._colortheme["COLOR_CAR_EGO"]
        else
            color = weighted_color_mean(1 - models[veh.id].c, colorant"red", colorant"green")
            # color_dict[veh.id] = RGB(HSV(0, 1.0 - models[veh.id].c, 1.0))
            color_dict[veh.id] = color
        # elseif models[veh.id].c == 1 && models[veh.id].idm.v_des < 15.0
        #     color_dict[veh.id] = cooperative #cooperative_slow
        # elseif models[veh.id].c == 0 && models[veh.id].idm.v_des < 15.0
        #     color_dict[veh.id] = aggressive # aggressive_slow
        # elseif models[veh.id].c == 1
        #     color_dict[veh.id] = cooperative
        # elseif models[veh.id].c == 0
        #     color_dict[veh.id] = aggressive
        end

    end
    return color_dict
end

struct BeliefOverlay <: SceneOverlay
    b::MergingBelief
end

function AutoViz.render!(rendermodel::RenderModel, overlay::BeliefOverlay,
                        scene::Scene, roadway::Roadway)
    for (vehid,prob) in overlay.b.driver_types
        vehind = findfirst(vehid, scene)
        if vehind != nothing
            render_proba!(rendermodel, scene, roadway, vehind, prob)
        end
    end
end

function render_proba!(rendermodel::RenderModel, scene::Scene, roadway::Roadway, vehind::Int64, proba::Float64)
    veh = scene[vehind]
    A = get_front(veh)
    disp_point = A + VecSE2(-5.0, 4.0, 0.0)
    text_overlay = TextOverlay(text=[@sprintf("%1.2f", proba)],
                                pos=disp_point,
                                line_spacing = 0.15,
                                font_size = 17,
                                incameraframe=true)
    render!(rendermodel, text_overlay, scene, roadway)
end


text:
include("AV.jl")
include("ecoli.jl")
include("manipulations.jl")


text:
mutable struct TranslatedShape{K,S<:AbstractShape{K}} <: AbstractShape{K}
    s::S
    ∆::SFloat{K}
end

translate(s::AbstractShape{K}, ∆::AbsVecReal) where {K} = translate(s, SFloat{K}(∆))
translate(s::AbstractShape{K}, ∆::SFloat{K}) where {K} = TranslatedShape(s, ∆)
translate(s::TranslatedShape{K}, ∆::SFloat{K}) where {K} = TranslatedShape(s, s.∆+∆)

level(x::SFloat{K}, s::TranslatedShape{K}) where {K} = level(x-s.∆, s.s)
project(x::SFloat{K}, s::TranslatedShape{K}) where {K} = project(x-s.∆, s.s)
ndir(x::SFloat{K}, s::TranslatedShape{K}) where {K} = ndir(x-s.∆, s.s) # outward direction even for x inside s

max_stars_repo_name: wsshin/LevelSetShapes

max_issues_repo_name: wsshin/LevelSetShapes

max_forks_repo_name: wsshin/LevelSetShapes


text:
include("cv.jl")

# glass

# Source: UCI / Glass Identification
# # of classes: 6
# # of data: 214
# # of features: 9
# Files:
# glass.scale
# > C=32, bias=1
# vote: 0.6261682242990654
# ctmc: 0.6355140186915887
# mean: 0.6355140186915887
# prod: 0.6261682242990654
features, labels = loaddata("data/glass.scale", 9)
cv(features, labels, C=2^5, bias=1)

# iris

# Source: UCI / Iris Plant
# # of classes: 3
# # of data: 150
# # of features: 4
# Files:
# iris.scale
# > C=32, bias=1
# vote: 0.98
# ctmc: 0.98
# mean: 0.98
# prod: 0.98
features, labels = loaddata("data/iris.scale", 4)
cv(features, labels, C=2^5, bias=1)

# segment

# Source: Statlog / Segment
# # of classes: 7
# # of data: 2,310
# # of features: 19
# Files:
# segment.scale
# > C=16, bias=1
# vote: 0.9575757575757575
# ctmc: 0.9575757575757575
# mean: 0.9571428571428572
# prod: 0.9571428571428572
features, labels = loaddata("data/segment.scale", 19)
cv(features, labels, C=2^5, bias=1)

# vehicle

# Source: Statlog / Vehicle
# # of classes: 4
# # of data: 846
# # of features: 18
# Files:
# vehicle.scale
# > C=32, bias=1
# vote: 0.8226950354609929
# ctmc: 0.8250591016548463
# mean: 0.8191489361702128
# prod: 0.8191489361702128
features, labels = loaddata("data/vehicle.scale", 18)
cv(features, labels, C=2^5, bias=1)

# vowel

# Source: UCI / Vowel
# Preprocessing: First 528 instances are used as training and the remaining instances are for testing. Scaling training data first and adjust testing data accordingly.
# # of classes: 11
# # of data: 528 / 462 (testing)
# # of features: 10
# Files:
# vowel
# vowel.t (testing)
# vowel.scale (scaled to [-1,1])
# vowel.scale.t (testing)
# > C=8, bias=1
# vote: 0.7537878787878788
# ctmc: 0.7424242424242424
# mean: 0.7443181818181818
# prod: 0.7443181818181818
features, labels = loaddata("data/vowel.scale", 10)
cv(features, labels, C=2^5, bias=1)

# wine

# Source: UCI / Wine Recognition
# # of classes: 3
# # of data: 178
# # of features: 13
# Files:
# wine.scale
# > C=0.0625, bias=1
# vote: 0.9775280898876404
# ctmc: 0.9775280898876404
# mean: 0.9775280898876404
# prod: 0.9775280898876404
features, labels = loaddata("data/wine.scale", 13)
cv(features, labels, C=2.0^-2, bias=1)

#####################

# dna

# Source: Statlog / Dna
# Preprocessing: Training data is further separated into two sets, tr and val. [CWH01a]
# # of classes: 3
# # of data: 2,000 / 1,186 (testing) / 1,400 (tr) / 600 (val)
# # of features: 180
# Files:
# dna.scale
# dna.scale.t (testing)
# dna.scale.tr (tr)
# dna.scale.val (val)
# > C=0.0078125, bias=4
# vote: 0.952
# ctmc: 0.9525
# mean: 0.9525
# prod: 0.9525
features, labels = loaddata("data/dna.scale", 180)
cv(features, labels, C=2.0^-7, bias=2^2)

# satimage

# Source: Statlog / Satimage
# Preprocessing: Training data is further separated into two sets, tr and val. [CWH01a]
# # of classes: 6
# # of data: 4,435 / 2,000 (testing) / 3,104 (tr) / 1,331 (val)
# # of features: 36
# Files:
# satimage.scale
# satimage.scale.t (testing)
# satimage.scale.tr (tr)
# satimage.scale.val (val)
# > C=8, bias=0.001
# vote: 0.8683201803833145
# ctmc: 0.86967305524239
# mean: 0.8683201803833145
# prod: 0.8683201803833145
features, labels = loaddata("data/satimage.scale", 36)
cv(features, labels, C=2^3, bias=.001)

# letter

# Source: Statlog / Letter
# Preprocessing: Training data is further separated into two sets, tr and val. [CWH01a]
# # of classes: 26
# # of data: 15,000 / 5,000 (testing) / 10,500 (tr) / 4,500 (val)
# # of features: 16
# Files:
# letter.scale
# letter.scale.t (testing)
# letter.scale.tr (tr)
# letter.scale.val (val)
# > C=8, bias=1
# vote: 0.8450666666666666
# ctmc: 0.8438666666666667
# mean: 0.8445333333333334
# prod: 0.8445333333333334
features, labels = loaddata("data/letter.scale", 16)
cv(features, labels, C=2^3, bias=1)

# shuttle

# Source: Statlog / Shuttle
# Preprocessing: Training data is further separated into two sets, tr and val. [CWH01a]
# # of classes: 7
# # of data: 43,500 / 14,500 (testing) / 30,450 (tr) / 13,050 (val)
# # of features: 9
# Files:
# shuttle.scale
# shuttle.scale.t (testing)
# shuttle.scale.tr (tr)
# shuttle.scale.val (val)
# > C=64, bias=1
# vote: 0.9693103448275862
# ctmc: 0.9695172413793104
# mean: 0.9717701149425287
# prod: 0.9717931034482759
features, labels = loaddata("data/shuttle.scale", 9)
cv(features, labels, C=2^6, bias=1)


text:
mutable struct SROCK1ConstantCache{zType,tType} <: StochasticDiffEqConstantCache
  ms::SVector{10,Int}
  mη::SVector{10,tType}
  zprev::zType
  mdeg::Int
  optimal_η::tType
end
@cache struct SROCK1Cache{uType,rateType,noiseRateType,C<:SROCK1ConstantCache} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  gₘ₋₁::noiseRateType
  gₘ₋₂::noiseRateType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  atmp::rateType
  constantcache::C
end

function alg_cache(alg::SROCK1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SROCK1ConstantCache{typeof(t)}(u)
end

function alg_cache(alg::SROCK1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  uᵢ₋₁ = zero(u)
  uᵢ₋₂ = zero(u)
  gₘ₋₁ = zero(noise_rate_prototype)
  gₘ₋₂ = zero(noise_rate_prototype)
  tmp  = uᵢ₋₂             # these 3 variables are dummied to use same memory
  fsalfirst = k
  atmp = zero(rate_prototype)
  constantcache = SROCK1ConstantCache{typeof(t)}(u)
  SROCK1Cache(u,uprev,uᵢ₋₁,uᵢ₋₂,gₘ₋₁,gₘ₋₂,tmp,k,fsalfirst,atmp,constantcache)
end

mutable struct SROCK2ConstantCache{zType,T} <: StochasticDiffEqConstantCache
  ms::SVector{46,Int}
  recf::Vector{T}
  mσ::SVector{46,T}
  mτ::SVector{46,T}
  recf2::Vector{T}
  mα::SVector{46,T}
  zprev::zType
  mdeg::Int
  deg_index::Int
  start::Int
end

@cache struct SROCK2Cache{uType,rateType,noiseRateType,T,C<:SROCK2ConstantCache} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ::uType
  uₓ::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  Gₛ::noiseRateType
  Gₛ₁::noiseRateType
  WikRange::T
  vec_χ::T
  tmp::uType
  k::rateType
  fsalfirst::rateType
  atmp::rateType
  constantcache::C
end

function alg_cache(alg::SROCK2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SROCK2ConstantCache{typeof(t)}(u)
end

function alg_cache(alg::SROCK2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  uᵢ = zero(u)
  uₓ = zero(u)
  uᵢ₋₁ = zero(u)
  uᵢ₋₂ = zero(u)
  Gₛ = zero(noise_rate_prototype)
  Gₛ₁ = zero(noise_rate_prototype)
  WikRange = false .* vec(ΔW)
  vec_χ = false .* vec(ΔW)
  tmp  = uᵢ₋₂            # these 2 variables are dummied to use same memory
  fsalfirst = k
  atmp = zero(rate_prototype)
  constantcache = SROCK2ConstantCache{typeof(t)}(u)
  SROCK2Cache(u,uprev,uᵢ,uₓ,uᵢ₋₁,uᵢ₋₂,Gₛ,Gₛ₁,WikRange,vec_χ,tmp,k,fsalfirst,atmp,constantcache)
end

mutable struct SROCKEMConstantCache{zType,tType} <: StochasticDiffEqConstantCache
  ms::SVector{10,Int}
  mη::SVector{10,tType}
  zprev::zType
  mdeg::Int
  optimal_η::tType
end

@cache struct SROCKEMCache{uType,rateType,noiseRateType,T,C<:SROCKEMConstantCache} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  Gₛ::noiseRateType
  Gₛ₁::noiseRateType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  WikRange::T
  atmp::rateType
  constantcache::C
end

function alg_cache(alg::SROCKEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SROCKEMConstantCache{typeof(t)}(u)
end

function alg_cache(alg::SROCKEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  uᵢ₋₁ = zero(u)
  uᵢ₋₂ = zero(u)
  Gₛ = zero(noise_rate_prototype)
  if (!alg.strong_order_1 || is_diagonal_noise(prob) || typeof(ΔW) <: Number || length(ΔW) == 1)
    Gₛ₁ = Gₛ
  else
    Gₛ₁ = zero(noise_rate_prototype)
  end
  WikRange = false .* vec(ΔW)
  tmp  = zero(u)             # these 3 variables are dummied to use same memory
  fsalfirst = k
  atmp = zero(rate_prototype)
  constantcache = SROCKEMConstantCache{typeof(t)}(u)
  SROCKEMCache(u,uprev,uᵢ₋₁,uᵢ₋₂,Gₛ,Gₛ₁,tmp,k,fsalfirst,WikRange,atmp,constantcache)
end

mutable struct SKSROCKConstantCache{zType,T} <: StochasticDiffEqConstantCache
  mc::Vector{T}
  mα::Vector{T}
  zprev::zType
end

@cache struct SKSROCKCache{uType,rateType,noise_rate_prototype,T,C<:SKSROCKConstantCache} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  Gₛ::noise_rate_prototype
  tmp::uType
  k::rateType
  fsalfirst::rateType
  WikRange::T
  atmp::rateType
  constantcache::C
end

function alg_cache(alg::SKSROCK,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SKSROCKConstantCache{typeof(t)}(u)
end

function alg_cache(alg::SKSROCK,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  uᵢ₋₁ = zero(u)
  uᵢ₋₂ = zero(u)
  Gₛ = zero(noise_rate_prototype)
  tmp  = uᵢ₋₂             # Dummmy variables
  fsalfirst = k
  WikRange = false .* vec(ΔW)
  atmp = zero(rate_prototype)
  constantcache = SKSROCKConstantCache{typeof(t)}(u)
  SKSROCKCache(u,uprev,uᵢ₋₁,uᵢ₋₂,Gₛ,tmp,k,fsalfirst,WikRange,atmp,constantcache)
end

mutable struct TangXiaoSROCK2ConstantCache{zType,T} <: StochasticDiffEqConstantCache
  ms::SVector{46,Int}
  recf::Vector{T}
  mσ::SVector{46,T}
  mτ::SVector{46,T}
  recf2::Vector{T}
  mα::SVector{5, T}
  mn̂::SVector{5, Int}
  c1::SVector{13, T}
  c2::SVector{13, T}
  zprev::zType
  mdeg::Int
  deg_index::Int
  start::Int
  start_mcs::Int
end

@cache struct TangXiaoSROCK2Cache{uType,rateType,noiseRateType,C<:TangXiaoSROCK2ConstantCache} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ::uType
  uₓ::uType
  Û₁::uType
  Û₂::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  Gₛ::noiseRateType
  Gₛ₁::noiseRateType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  atmp::rateType
  constantcache::C
end

function alg_cache(alg::TangXiaoSROCK2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  TangXiaoSROCK2ConstantCache{typeof(t)}(u)
end

function alg_cache(alg::TangXiaoSROCK2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  uᵢ = zero(u)
  uₓ = zero(u)
  Û₁ = zero(u)
  Û₂ = zero(u)
  uᵢ₋₁ = zero(u)
  uᵢ₋₂ = zero(u)
  Gₛ = zero(noise_rate_prototype)
  if typeof(ΔW) <: Number || length(ΔW) == 1 || is_diagonal_noise(prob)
    Gₛ₁ = Gₛ
  else
    Gₛ₁ = zero(noise_rate_prototype)
  end
  tmp  = uᵢ₋₂            # these 2 variables are dummied to use same memory
  fsalfirst = k
  atmp = zero(rate_prototype)
  constantcache = TangXiaoSROCK2ConstantCache{typeof(t)}(u)
  TangXiaoSROCK2Cache(u,uprev,uᵢ,uₓ,Û₁,Û₂,uᵢ₋₁,uᵢ₋₂,Gₛ,Gₛ₁,tmp,k,fsalfirst,atmp,constantcache)
end

mutable struct KomBurSROCK2ConstantCache{zType,T} <: StochasticDiffEqConstantCache
  ms::SVector{46,Int}
  recf::Vector{T}
  mσ::SVector{46,T}
  mτ::SVector{46,T}
  mδ::Vector{T}
  zprev::zType
  mdeg::Int
  deg_index::Int
  start::Int
end

@cache struct KomBurSROCK2Cache{uType,rateType,noiseRateType,T,C<:KomBurSROCK2ConstantCache} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  utmp::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  k::rateType
  yₛ₋₁::rateType
  yₛ₋₂::rateType
  yₛ₋₃::rateType
  SXₛ₋₁::uType
  SXₛ₋₂::uType
  SXₛ₋₃::uType
  Gₛ::noiseRateType
  Xₛ₋₁::noiseRateType
  Xₛ₋₂::noiseRateType
  Xₛ₋₃::noiseRateType
  vec_χ::T
  tmp::uType
  fsalfirst::rateType
  WikRange::T
  atmp::rateType
  constantcache::C
end

function alg_cache(alg::KomBurSROCK2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  KomBurSROCK2ConstantCache{typeof(t)}(u)
end

function alg_cache(alg::KomBurSROCK2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  utmp = zero(u)
  uᵢ₋₁ = zero(u)
  uᵢ₋₂ = zero(u)
  k = zero(rate_prototype)
  yₛ₋₁ = k
  yₛ₋₂ = zero(rate_prototype)
  yₛ₋₃ = zero(rate_prototype)
  Xₛ₋₁ = zero(noise_rate_prototype)
  Xₛ₋₂ = zero(noise_rate_prototype)
  Xₛ₋₃ = zero(noise_rate_prototype)
  vec_χ = false .* vec(ΔW)
  WikRange = false .* vec(ΔW)
  if typeof(ΔW) <: Number || length(ΔW) == 1 || is_diagonal_noise(prob)
    Gₛ = Xₛ₋₁
    SXₛ₋₁ = utmp
    SXₛ₋₂ = utmp
    SXₛ₋₃ = utmp
  else
    Gₛ = zero(noise_rate_prototype)
    SXₛ₋₁ = zero(uprev)
    SXₛ₋₂ = zero(uprev)
    SXₛ₋₃ = zero(uprev)
  end
  tmp  = uᵢ₋₂            # these 3 variables are dummied to use same memory
  fsalfirst = k
  atmp = yₛ₋₂
  constantcache = KomBurSROCK2ConstantCache{typeof(t)}(u)
  KomBurSROCK2Cache(u,uprev,utmp,uᵢ₋₁,uᵢ₋₂,k,yₛ₋₁,yₛ₋₂,
                    yₛ₋₃,SXₛ₋₁,SXₛ₋₂,SXₛ₋₃,Gₛ,Xₛ₋₁,Xₛ₋₂,Xₛ₋₃,vec_χ,
                    tmp,fsalfirst,WikRange,atmp,constantcache)
end

mutable struct SROCKC2ConstantCache{zType,T} <: StochasticDiffEqConstantCache
  ms::SVector{46,Int}
  recf::Vector{T}
  mσ::SVector{46,T}
  mτ::SVector{46,T}
  zprev::zType
  mdeg::Int
  deg_index::Int
  start::Int
end

@cache struct SROCKC2Cache{uType,rateType,noiseRateType,T,C<:SROCKC2ConstantCache} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uᵢ::uType
  uᵢ₋₁::uType
  uᵢ₋₂::uType
  Gₛ::noiseRateType
  Gₛ₁::noiseRateType
  WikRange::T
  tmp::uType
  k::rateType
  fsalfirst::rateType
  atmp::rateType
  constantcache::C
end

function alg_cache(alg::SROCKC2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SROCKC2ConstantCache{typeof(t)}(u)
end

function alg_cache(alg::SROCKC2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  k = zero(rate_prototype)
  uᵢ = zero(u)
  uᵢ₋₁ = zero(u)
  uᵢ₋₂ = zero(u)
  Gₛ = zero(noise_rate_prototype)
  Gₛ₁ = zero(noise_rate_prototype)
  WikRange = false .* vec(ΔW)
  tmp  = zero(u)
  fsalfirst = k           # this variables are dummied to use same memory
  atmp = zero(rate_prototype)
  constantcache = SROCKC2ConstantCache{typeof(t)}(u)
  SROCKC2Cache(u,uprev,uᵢ,uᵢ₋₁,uᵢ₋₂,Gₛ,Gₛ₁,WikRange,tmp,k,fsalfirst,atmp,constantcache)
end


text:
@testset "moduleprovider" begin

Context() do ctx
let
    mod = LLVM.Module("SomeModule"; ctx)
    mp = ModuleProvider(mod)
    dispose(mp)
end
end

Context() do ctx
let
    mod = LLVM.Module("SomeModule"; ctx)
    ModuleProvider(mod) do mp

    end
end
end

end


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

"""
    MulticlassDistance()(x, y)

Return `true` if x and y are different classes, else return `false`.
"""
struct MulticlassDistance <: Metric end

(::MulticlassDistance)(x, y) = x != y

result_type(::MulticlassDistance, x, y) = Bool

"""
    OrderedFactorDistance()(x, y)

Return the absolute value of the difference between the categorical codes of x and y.
"""
struct OrderedFactorDistance <: Metric end

(::OrderedFactorDistance)(x, y) = abs(levelcode(x) - levelcode(y))

result_type(::OrderedFactorDistance, x, y) = Float64

default_distance(::Type{Continuous})      = Distances.Euclidean()
default_distance(::Type{Count})           = Distances.Cityblock()
default_distance(::Type{<:Multiclass})    = MulticlassDistance()
default_distance(::Type{<:OrderedFactor}) = OrderedFactorDistance()
default_distance(::Type{Textual})         = StringDistances.Levenshtein()
default_distance(::Type{<:Compositional}) = CoDa.Aitchison()

function default_distances(table)
  columns = Tables.columnnames(table)
  scitypes = schema(table).scitypes
  distances = default_distance.(scitypes)
  Dict(columns .=> distances)
end


text:
function __init__()
    RLDatasets.d4rl_init()
    RLDatasets.d4rl_pybullet_init()
    RLDatasets.atari_init()
    RLDatasets.rl_unplugged_atari_init()
end