text:
struct _Undefined end
struct _Defined end

"""
    rasterize(x; kw...)
    rasterize(points, values; kw...)

Rasterize the points and values in a Tables.jl compatible object,
or separate `points` and `values` vetors/iterators, which must be
the same length.

If a GeoInterface `AbstractGeometry` or nested `Vector`s of `Tuple`/`Vector`
points is passed in, a `fill` keyword is also required to provide the value
that the array will be filled with. If `fill` is a function, it will be
applied to the existing value present in the array.

# Arguments

- `x`: a Tables.jl compatible object containing points and values columns,
    an GeoInterface.jl `AbstractGeometry`, or a nested `Vector` of `Vectors`.
- `points`: A `Vector` or nested `Vectors` holding `Vector` or `Tuple` of `Real`
- `values` A `Vector` of values to be written to a `Raster`, or a Vector of `NamedTuple`
    to write to a `RasterStack`.

# Keywords

These are detected automatically from `A` and `data` where possible.

- `to`: a `Raster`, `RasterStack` of `Tuple` of `Dimension` to use as a to.
- `order`: A `Tuple` of pairs `Dim => Symbol` for the keys in the data that match
    the dimension, or for the order of dimensions in ppoint values, like `(X, Y)`.
- `atol`: an absolute tolerance for rasterizing to dimensions with `Points` sampling.
- `filename`: a filename to write to directly, useful for large files.
- `suffix`: a string or value to append to the filename.
    A tuple of `suffix` will be applied to stack layers. `keys(st)` are the default.

## Geometry keywords

These can be used when a `GeoInterface.AbstractGeometry` is passed in.

- `fill`: the value to fill a polygon with, if `data` is a polygon.
- `shape`: Force `data` to be treated as `:polygon`, `:line` or `:point`.

And specifically for `shape=:polygon`:

- `boundary`: include pixels where the `:center` is inside the polygon, where
    the line `:touches` the pixel, or that are completely `:inside` inside the polygon.

## Table keywords

- `name`: A `Symbol` to return a `Raster` from a single column,
    or `Tuple` of `Symbol` to return a `RasterStack` from multiple columns.

# Example

Rasterize a shapefile for China and plot, with a border.

```jldoctest
using Rasters, Plots, Dates, Shapefile, Downloads
using Rasters.LookupArrays

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary_lines.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)

# Loade the shapes for china
china_border = Shapefile.Handle(shapefile_name).shapes[10]

# Define dims for the china area
dms = Y(Projected(15.0:0.1:55.0; order=ForwardOrdered(), span=Regular(0.1), sampling=Intervals(Start()), crs=EPSG(4326))),
      X(Projected(70.0:0.1:140; order=ForwardOrdered(), span=Regular(0.1), sampling=Intervals(Start()), crs=EPSG(4326)))

# Rasterize the border polygon
china = rasterize(china_border;
    to=dms, missingval=-9999, fill=1,
    order=(X, Y), shape=:polygon,
    boundary=:touches,
)

# And plot
p = plot(china; color=:spring)
plot!(p, china_border; fillalpha=0, linewidth=0.6)
savefig("build/china_rasterized.png")

# output
```

![rasterize](china_rasterized.pngfill)

$EXPERIMENTAL
"""
function rasterize(data; to, order=nothing, name=nothing, kw...)
    if Tables.istable(data)
        points, values, order, name = _table_to_points_values(to, data; order, name)
        return _rasterize(to, points, values; order, name, kw...)
    else
        order = isnothing(order) ? DEFAULT_POINT_ORDER : order
        order = dims(to, order)
        return _rasterize(to, data; order, init=_Undefined(), kw...)
    end
end
rasterize(points, values; to, kw...)  = _rasterize(to, points, values; kw...)

function _rasterize(to::DimTuple, points;
    filename=nothing, suffix=nothing, metadata=NoMetadata(), keys=nothing, name=keys, parent=nothing,
    fill, eltype=typeof(fill), missingval=_writeable_missing(filename, eltype), kw...
)
    A = _alloc_rasterize(filename, eltype, to; missingval, suffix, parent) do a
        rasterize!(a, points; fill, missingval, name, kw...)
    end
    return A
end
function _rasterize(to::DimTuple, points, vals;
    missingval=nothing, filename=nothing, suffix=nothing,
    metadata=NoMetadata(), keys=nothing, name=keys, parent=nothing, kw...
)
    firstval = first(vals)

    if firstval isa Union{NamedTuple,Tuple,AbstractArray}
        keys = _filter_name(name, vals)
        # Rasterize mutiple values to a stack
        layers = map(keys, values(firstval)) do key, val
            missingval = isnothing(missingval) ? _writeable_missing(filename, typeof(val)) : missingval
            _alloc_rasterize(filename, typeof(val), to; name, metadata, missingval, suffix=key, parent) do a
                a .= missingval
            end
        end |> NamedTuple{keys}
        st = RasterStack(layers, to)
        return rasterize!(st, points, vals; kw...)
    else
        name = first(cleankeys(name))
        # Rasterize to an array
        missingval = isnothing(missingval) ? _writeable_missing(filename, typeof(firstval)) : missingval
        A = _alloc_rasterize(filename, typeof(firstval), to; name, metadata, missingval, suffix, parent) do a
            a .= missingval
        end
        return rasterize!(A, points, vals; kw...)
    end
end
function _rasterize(to::AbstractRaster, args...;
    missingval=missingval(to), keys=name(to), name=keys, kw...
)
    _rasterize(dims(to), args...; name, missingval, kw...)
end
function _rasterize(to::AbstractRasterStack, args...; keys=keys(to), name=keys, kw...)
    return _rasterize(dims(to), args...; name, kw...)
end

function _alloc_rasterize(f, filename, T, to; missingval, suffix=nothing, kw...)
    T = promote_type(typeof(missingval), T)
    A = create(filename, T, to; suffix, missingval, kw...)
    open(f, A; write=true)
    return A
end

"""
    rasterize!(x, data; order, name, atol)
    rasterize!(x, points, vals; order, atol)

Rasterize the points and vals in `data`, or the `points` and `vals` objects,
into the [`Raster`](@ref) or [`RasterStack`](@ref) `x`.

# Arguments

- `x`: a `Raster` or `RasterStack` to rasterize to.
- `data`: a Tables.jl compatible object containing points and values or a
    polygon - an GeoInterface.jl `AbstractGeometry`, or a nested `Vector` of `Vectors`.
- `points`: A `Vector` or nested `Vector` holding `Vector` or `Tuple` of `Real`
- `vals` A `Vector` of values to be written when `x` is a `Raster`, or a Vector of
    `NamedTuple` to write when `x` is a `RasterStack`.

# Keywords

These are detected automatically from `A` and `data` where possible.

- `order`: A `Tuple` of pairs `Dim => Symbol` for the keys in the data that match
    the dimension, or for the order of dimensions in points, like `(X, Y)`.
- `atol`: an absolute tolerance for rasterizing to dimensions with `Points` sampling.
- `filename`: a filename to write to directly, useful for large files.
- `suffix`: a string or value to append to the filename.
    A tuple of `suffix` will be applied to stack layers. `keys(st)` are the default.

## Geometry keywords

These can be used when a `GeoInterface.AbstractGeometry` is passed in.

- `fill`: the value to fill a polygon with, if `data` is a polygon.
    `fill can also be a `Function` of the existing value.
- `shape`: Force `data` to be treated as `:polygon`, `:line` or `:point`.

And specifically for `shape=:polygon`:

- `boundary`: include pixels where the `:center` is inside the polygon, where
    the line `:touches` the pixel, or that are completely `:inside` inside the polygon.

## Table keywords

- `name`: A `Symbol` to return a `Raster` from a single column,
    or `Tuple` of `Symbol` to return a `RasterStack` from multiple columns.

# Example

```jldoctest
using Rasters, Plots, Dates, Shapefile, GeoInterface, Downloads
using Rasters.LookupArrays

# Download a borders shapefile
shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "indonesia_border.shp"
isfile(shapefile_name) || Downloads.download(shapefile_url, shapefile_name)

# Load the shapes for denmark
indonesia_border = Shapefile.Handle(shapefile_name).shapes[1]

# Make an empty EPSG 4326 projected Raster of the area of Indonesia
dimz = Y(-15.0:0.1:10.9; mode=Projected(; sampling=Intervals(Start()), crs=EPSG(4326))),
       X(90.0:0.1:145; mode=Projected(; sampling=Intervals(Start()), crs=EPSG(4326)))

A = Raster(zeros(UInt16, dimz); missingval=0)

# Rasterize each island with a different number
for (i, shp) in enumerate(coordinates(indonesia_border))
    rasterize!(A, shp; fill=i, order=(X, Y), shape=:polygon, boundary=:touches)
end

# And plot
p = plot(A; color=:spring)
plot!(p, indonesia_border; fillalpha=0, linewidth=0.7)
savefig("build/indonesia_rasterized.png")

# output

```

![rasterize](indonesia_rasterized.png)

$EXPERIMENTAL
"""
function rasterize!(x::RasterStackOrArray, data;
    order=nothing, name=nothing, missingval=missingval(x), kw...
)
    if Tables.istable(data)
        _rasterize_table!(x, data; name, kw...)
    else
        order = isnothing(order) ? DEFAULT_POINT_ORDER : order
        buffer = Raster(falses(commondims(x, order)))
        _rasterize_geometry!(x, data; order, buffer, init=_Defined(), missingval, kw...)
    end
end
function rasterize!(A::AbstractRaster, geom::GI.AbstractGeometry, vals; shape=:point, kw...)
    rasterize!(A, _flat_nodes(geom), vals; kw...)
end
function rasterize!(st::AbstractRasterStack, geom::GI.AbstractGeometry, vals; shape=:point, kw...)
    rasterize!(st, _flat_nodes(geom), vals; kw...)
end
function rasterize!(st::AbstractRasterStack, points, vals; keys=keys(st), name=keys, kw...)
    keys = _filter_name(name, vals)
    return _rasterize!(st[name], points, vals; name, kw...)
end
function rasterize!(A::AbstractRaster, points, vals; name=nothing, kw...)
    _rasterize!(A, points, vals; name, kw...)
end

function _rasterize!(x::RasterStackOrArray, points, vals;
    order=DEFAULT_POINT_ORDER, atol=nothing, name=nothing
)
    check_points(points, order)
    isdisk(first(x)) && _warn_disk()
    ordered_dims = dims(x, order)
    _without_mapped_crs(x) do x
        map(points, vals) do p, v
            any(map(ismissing, p)) && return nothing
            selectors = map((d, x) -> _at_or_contains(d, x, atol), ordered_dims, p)
            all(map(DD.hasselection, ordered_dims, selectors)) || return nothing
            _set_at_selection!(x, selectors, v)
            return nothing
        end
    end
    return x
end

function _set_at_selection!(st::AbstractRasterStack, selectors, vals)
    map(values(st), values(vals)) do A, val
        _set_at_selection!(A, selectors, val)
    end
end
function _set_at_selection!(A::AbstractRaster, selectors, val::Union{NamedTuple,Tuple})
    _set_at_selection!(A, selectors, first(val))
end
function _set_at_selection!(A::AbstractRaster, selectors, val)
    if length(selectors) == length(dims(A))
        A[selectors...] = val
    else
        A[selectors...] .= val
    end
end

function check_points(points, order)
    lp = length(first(points))
    lo = length(order)
    if lp !== lo
        throw(ArgumentError("length of points in `points` ($lp) is different to the length of `order` ($lo)"))
    end
end

# _rasterize_table!
function _rasterize_table!(x::RasterStackOrArray, data; order=nothing, name=nothing, kw...)
    points, vals, order, name = _table_to_points_values(x, data; order, name)
    return rasterize!(x, points, vals; order, name, kw...)
end

# _rasterize_geometry!
function _rasterize_geometry!(x::RasterStackOrArray, geom;
    buffer, fill, order, init, missingval, kw...
)
    boolmask!(buffer, geom; order, kw...)
    return _fill!(x, buffer, fill, init, missingval)
end

_fill(B, fill, missingval) where T = broadcast_dims(x -> x ? fill : missingval, B)

function _fill!(st::AbstractRasterStack, B, fill, args...)
    map((a, f) -> _fill!(a, B, f, args...), st, fill)
    return st
end
# If the array is initialised, we can use the existing values
function _fill!(A::AbstractRaster{T}, B, fill, init::_Defined, missingval) where T
    broadcast_dims!(A, A, B) do a, b
        val = b ? (fill isa Function ? fill(a) : fill) : a
        convert(T, val) # In case we are writing to disk
    end
end
# If the array is not yet initialised, we have to fill with fill and misssingval
function _fill!(A::AbstractRaster{T}, B, fill, init::_Undefined, missingval) where T
    fill = convert(T, fill) # In case we are writing to disk
    missingval = convert(T, missingval)
    broadcast_dims!(A, B) do b
        b ? fill : missingval
    end
end

function _at_or_contains(d, v, atol)
    selector = sampling(d) isa Intervals ? Contains(v) : At(v; atol=atol)
    DD.basetypeof(d)(selector)
end

function _filter_name(name, vals)
    firstval = first(vals)
    name = if isnothing(name)
        if firstval isa NamedTuple
            keys(firstval)
        else
            stackkeys
        end
    else
        if firstval isa NamedTuple
            name == keys(firstval) || throw(ArgumentError("`name` keyword does not match point names"))
        elseif firstval isa Union{Tuple,Array}
            length(name) == length(firstval) || throw(ArgumentError("`name` keyword does not match length of point"))
        end
        name
    end
end


text:
@testset "Structure tests" begin

    #Helper functions
    nt(name) = t -> NamedTrait(t, name)
    tf(name) = t -> nt(name)(t)
    inpt(size, id="in") = inputvertex(id, size)
    av(in, outsize; name="av", comp = identity) = absorbvertex(comp, outsize, in, traitdecoration = tf(name))
    sv(in, ins...; name="sv") = conc(in, ins..., dims=2, traitdecoration = tf(name))
    sv(in; name="sv") = vertex(identity, nout(in), tf(name)(SizeStack()), in)
    iv(ins...; name="iv") = +(traitconf(tf(name)) >> ins[1], ins[2:end]...)
    imu(in, outsize; name="imu") = immutablevertex(identity, outsize, in, traitdecoration= tf(name))

    @testset "Edge mutation" begin
        @testset "Edge removal" begin

            @testset "Fail noop" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v1, 4, name="v2")

                @test remove_edge!(v1, v2; strategy=FailAlignSizeNoOp()) == false

                @test name.(inputs(v2)) == [name(v1)]
                @test name.(outputs(v1)) == [name(v2)]
            end

            @testset "Remove from absorbing" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v1, 7, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                #Only because helper method av does not take muliple inputs!
                create_edge!(v2, v3)
                @test inputs(v3) == [v1, v2]

                #Now for the actual test...
                remove_edge!(v2, v3)
                @test inputs(v3) == [v1]
                @test outputs(v1) == [v3]
                @test nin(v3) == [nout(v1)] == [5]

                @test inputs(v5) == [v2]
                @test outputs(v2) == [v5]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Remove from stacking" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = sv(v1, v2, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                remove_edge!(v2, v3)
                @test inputs(v3) == [v1]
                @test outputs(v1) == [v3]
                @test nin(v4) == [nout(v3)] == nin(v3) == [nout(v1)] == [5]

                @test inputs(v5) == [v2]
                @test outputs(v2) == [v5]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Remove from stacking and increase" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = sv(v1, v2, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                remove_edge!(v2, v3, strategy=PostAlignJuMP(ΔNoutExact(v3, 2)))
                @test inputs(v3) == [v1]
                @test outputs(v1) == [v3]
                @test nin(v4) == [nout(v3)] == nin(v3) == [nout(v1)] == [7]

                @test inputs(v5) == [v2]
                @test outputs(v2) == [v5]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Remove duplicate from single output stacking" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = sv(v1, v2, v1, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1, v2, v1]
                remove_edge!(v1, v3)

                @test inputs(v3) == [v2, v1]
                @test nin(v4) == [nout(v3)] == [nout(v1) + nout(v2)] == [9]
                @test nin(v3) == [nout(v2), nout(v1)] == [4, 5]

                @test outputs(v2) == [v3, v5]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Remove other duplicate from single output stacking" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = sv(v1, v2, v1, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1, v2, v1]
                remove_edge!(v1, v3, nr=2)

                @test inputs(v3) == [v1, v2]
                @test nin(v4) == [nout(v3)] == [nout(v1) + nout(v2)] == [9]
                @test nin(v3) == [nout(v1), nout(v2)] == [5, 4]

                @test outputs(v2) == [v3, v5]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Remove from invariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 8, name="v1")
                v2 = av(v0, 8, name="v2")
                v3 = iv(v1, v2, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                remove_edge!(v2, v3)
                @test inputs(v3) == [v1]
                @test outputs(v1) == [v3]
                @test nin(v3) == [nout(v1)] == [8]

                @test inputs(v5) == [v2]
                @test outputs(v2) == [v5]
                @test nin(v5) == [nout(v2)] == [8]
            end

            @testset "Remove duplicate from single output invariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = iv(v1, v2, v1, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1, v2, v1]
                remove_edge!(v1, v3)

                @test inputs(v3) == [v2, v1]
                @test nin(v4) == [nout(v3)] == [nout(v1)] == [nout(v2)] == [4]

                @test outputs(v2) == [v3, v5]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Remove other duplicate from single output invariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = iv(v1, v2, v1, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1, v2, v1]
                remove_edge!(v1, v3, nr=2)

                @test inputs(v3) == [v1, v2]
                @test nin(v4) == [nout(v3)] == [nout(v1)] == [nout(v2)] == [4]

                @test outputs(v2) == [v3, v5]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Revert remove edge SizeStack" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = sv(v0, v1, v2, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v0, v1, v2]
                @test nin(v3) == nout.([v0, v1, v2]) == [3,5,4]
                @test [nout(v3)] == nin(v4) == [3+5+4]
                @test nin_org(v3) == nout_org.([v0, v1, v2]) == [3,5,4]
                @test [nout_org(v3)] == nin_org(v4) == [3+5+4]

                struct RevertPost <: AbstractAlignSizeStrategy end
                NaiveNASlib.postalignsizes(::RevertPost, vin, vout, pos) = NaiveNASlib.postalignsizes(FailAlignSizeRevert(), vin, vout, pos)

                remove_edge!(v1, v3, strategy=RevertPost())

                @test inputs(v3) == [v0, v1, v2]
                @test nin(v3) == nout.([v0, v1, v2]) == [3,5,4]
                @test [nout(v3)] == nin(v4) == [3+5+4]
                @test nin_org(v3) == nout_org.([v0, v1, v2]) == [3,5,4]
                @test [nout_org(v3)] == nin_org(v4) == [3+5+4]
            end

            @testset "PostSelectOutputs SizeInvariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 4, name="v3")
                v4 = iv(v1, v2, v3, name = "v4")
                v5 = av(v4, 3, name="v5")

                remove_edge!(v3, v4, strategy=PostSelectOutputs(align = DecreaseBigger(), valuefun = v -> 1:nout_org(v)))

                @test inputs(v4) == [v1, v2]
                @test nin(v4) == nout.([v1, v2]) == [4,4]
                @test in_inds(op(v4)) == out_inds.(op.([v1,v2])) == [1:4, 1:4]
            end

            @testset "PostSelectOutputs SizeInvariant post align" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 4, name="v3")
                v4 = iv(v1,v2,v3, name = "v4")
                v5 = av(v4, 3, name="v5")

                remove_edge!(v3, v4, strategy=PostSelectOutputs(valuefun = v -> 1:nout_org(v)))

                @test inputs(v4) == [v1, v2]
                @test nin(v4) == nout.([v1, v2]) == [4,4]
                @test in_inds(op(v4)) == out_inds.(op.([v1,v2])) == [1:4, 1:4]
            end

            @testset "PostSelectOutputs SizeStack" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = sv(v1,v2,v3, name = "v4")
                v5 = av(v4, 3, name="v5")

                remove_edge!(v2, v4, strategy=PostSelectOutputs(valuefun = v -> 1:nout_org(v)))

                @test inputs(v4) == [v1, v3]
                @test nin(v4) == nout.([v1, v3]) == [3,5]
                @test out_inds(op(v1)) == 1:3
                @test out_inds(op(v2)) == 1:4
                @test out_inds(op(v3)) == 1:5
                # This would be better if it was [1:3;8:12] but remove_vertex removes the edge before PostSelectOutputs has a chance to see it :(
                @test out_inds(op(v4)) == 1:8
                @test in_inds(op(v5)) == [1:8]
            end

            @testset "PostSelectOutputs fail size align" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = sv(v1,v2,v3, name = "v4")
                v5 = av(v4, 3, name="v5")

                @test_logs (:warn, r"Could not align size") remove_edge!(v1, v4, strategy=PostSelectOutputs(align = FailAlignSizeWarn()))

                @test inputs(v4) == [v1, v2, v3]
                @test nin(v4) == nout.([v1, v2, v3]) == [3,4,5]
                @test nin_org(v4) == nout_org.([v1, v2, v3]) == [3,4,5]
                @test [nout(v4)] == nin(v5) == [3+4+5]
                @test [nout_org(v4)] == nin_org(v5) == [3+4+5]
            end

            @testset "PostSelectOutputs fail select" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = sv(v1,v2,v3, name = "v4")
                v5 = av(v4, 3, name="v5")

                @test_logs (:warn, r"Could not align size") remove_edge!(v2, v4, strategy=PostSelectOutputs(
                align = PostAlignJuMP(), select = NoutRevert(), fallback=FailAlignSizeWarn()))

                @test inputs(v4) == [v1, v2, v3]
                @test nin(v4) == nout.([v1, v2, v3]) == [3,4,5]
                @test nin_org(v4) == nout_org.([v1, v2, v3]) == [3,4,5]
                @test [nout(v4)] == nin(v5) == [3+4+5]
                @test [nout_org(v4)] == nin_org(v5) == [3+4+5]
            end

            @testset "PostApplyMutation SizeStack" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = sv(v1,v2,v3, name = "v4")
                v5 = av(v4, 3, name="v5")

                remove_edge!(v1, v4, strategy=PostApplyMutation())

                @test inputs(v4) == [v2, v3]
                @test nin(v4) == nin_org(v4) == nout.([v2, v3]) == [4,5]
                @test nin(v5) == nin_org(v5) == [nout(v4)] == [nout_org(v4)] == [4+5]
            end

            @testset "Remove all $t" for (t, vf) in ((SizeStack, sv), (SizeInvariant, iv))
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 3, name="v2")
                v3 = av(v0, 3, name="v3")
                v4 = vf(v1,v2,v3, name = "v4")
                v5 = av(v4, 3, name="v5")
                v6 = av(v4, 4, name="v6")

                is = copy(inputs(v4))
                os = copy(outputs(v4))

                @test all(ov -> remove_edge!(v4, ov; strategy=NoSizeChange()), os)
                @test all(iv -> remove_edge!(iv, v4; strategy=NoSizeChange()), is)

                @test inputs(v4) == []
                @test outputs(v4) == []

                @test outputs(v1) == []
                @test outputs(v2) == []
                @test outputs(v3) == []

                @test inputs(v5) == []
                @test inputs(v6) == []
            end
        end

        @testset "Edge addition" begin

            @testset "Fail noop" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v1, 4, name="v2")

                @test create_edge!(v1, v2; strategy=FailAlignSizeNoOp()) == false

                @test name.(inputs(v2))== [name(v1)]
                @test name.(outputs(v1)) == [name(v2)]
            end

            @testset "Add to absorbing" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v1, 7, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1]
                create_edge!(v2, v3)

                @test inputs(v3) == [v1, v2]
                @test nin(v3) == [nout(v1), nout(v2)] == [5, 4]
                @test nin(v4) == [nout(v3)] == [7]

                @test outputs(v2) == [v5, v3]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Add to single output stacking" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = sv(v1, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1]
                create_edge!(v2, v3)

                @test inputs(v3) == [v1, v2]
                @test nin(v4) == [nout(v3)] == [nout(v1) + nout(v2)] == [9]

                @test outputs(v2) == [v5, v3]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Add to single output stacking and increase" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = sv(v1, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1]
                create_edge!(v2, v3, strategy=PostAlignJuMP(ΔNoutExact(v2, 2)))

                @test inputs(v3) == [v1, v2]
                @test nin(v4) == [nout(v3)] == [nout(v1) + nout(v2)] == [11]

                @test outputs(v2) == [v5, v3]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [6]
            end

            @testset "Add duplicate to single output stacking" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = sv(v1, v1, v2, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1, v1, v2]
                create_edge!(v1, v3)

                @test inputs(v3) == [v1, v1, v2, v1]
                @test nin(v4) == [nout(v3)] == [3nout(v1) + nout(v2)] == [19]

                @test outputs(v2) == [v3, v5]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Add immutable to single output stacking" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = sv(v1, name = "v2")
                v3 = av(v2, 3, name="v3")

                @test inputs(v2) == [v1]
                create_edge!(v0, v2)

                @test inputs(v2) == [v1, v0]
                @test nin(v3) == [nout(v2)] == [nout(v1) + nout(v0)] == [8]

                @test outputs(v0) == [v1, v2]
            end

            @testset "Add to single output invariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = iv(v1, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1]
                create_edge!(v2, v3)

                @test inputs(v3) == [v1, v2]
                @test nin(v4) == [nout(v3)] == [nout(v1)] == [nout(v2)] == [5]

                @test outputs(v2) == [v5, v3]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [5]
            end

            @testset "Add duplicate to single output invariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = iv(v1, v1, v2, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = av(v2, 2, name="v5")

                @test inputs(v3) == [v1, v1, v2]
                create_edge!(v1, v3)

                @test inputs(v3) == [v1, v1, v2, v1]
                @test nin(v4) == [nout(v3)] == [nout(v1)] == [nout(v2)] == [4]

                @test outputs(v2) == [v3, v5]
                @test inputs(v5) == [v2]
                @test nin(v5) == [nout(v2)] == [4]
            end

            @testset "Add immutable to single output invariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = iv(v1, name = "v2")
                v3 = av(v2, 3, name="v3")

                @test inputs(v2) == [v1]
                create_edge!(v0, v2)

                @test inputs(v2) == [v1, v0]
                @test nin(v3) == [nout(v2)] == [nout(v1)] == [nout(v0)] == [3]

                @test outputs(v0) == [v1, v2]
            end

            @testset "Add with hidden SizeStack" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 4, name="v2")
                vh = av(v0, 5, name="vh")
                v3 = sv(v1, name = "v3")
                v4 = av(v3, 3, name="v4")
                v5 = sv(v4, vh, name="v5")
                v6 = av(v2, 2, name="v6")

                @test inputs(v3) == [v1]
                create_edge!(v2, v3)

                @test inputs(v3) == [v1, v2]
                @test nin(v4) == [nout(v3)] == [nout(v1) + nout(v2)] == [9]

                @test outputs(v2) == [v6, v3]
                @test inputs(v6) == [v2]
                @test nin(v6) == [nout(v2)] == [4]

                @test outputs(vh) == [v5]
                @test inputs(v5) == [v4, vh]
                @test nin(v5) == [nout(v4), nout(vh)] == [3, 5]
            end

            @testset "Add with hidden SizeInvariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 5, name="v1")
                v2 = av(v0, 3, name="v2")
                vh = av(v0, 5, name="vh")
                v3 = iv(v1, vh, name = "v3")
                v4 = av(v3, 3, name="v4")

                @test inputs(v3) == [v1, vh]
                create_edge!(v2, v3, strategy=AlignSizeBoth())

                @test inputs(v3) == [v1, vh, v2]
                @test nin(v4) == [nout(v3)] == [nout(v1)] == [nout(vh)] == [nout(v2)] == [4]
            end

            @testset "Add with size cycle SizeInvariant" begin
                v1 = av(inpt(3, "in"), 3, name="v1")
                v2 = av(v1, 5, name="v2")
                v3 = iv(v2, name="v3")
                p1 = av(v3, 3, name="p1")
                p2 = av(v3, 2, name="p2")
                v4 = sv(p1,p2, name="v4")
                v5 = "v5" >> v4 + v2

                @test inputs(v5) == [v4, v2]
                @test nin(v5) == [nout(v4), nout(v2)] == [5, 5]

                @test_logs (:warn, r"Can not add edge") create_edge!(p1, v5)
                @test inputs(v5) == [v4, v2]
                @test nin(v5) == [nout(v4), nout(v2)] == [5, 5]
            end

            @testset "Add with size cycle SizeStack" begin
                v1 = av(inpt(3, "in"), 3, name="v1")
                v2 = av(v1, 5, name="v2")
                v3 = iv(v2, name="v3")
                ve = iv(v1, name="ve")
                v4 = sv(v3, name="v4")
                v5 = "v5" >> v4 + v2

                @test inputs(v5) == [v4, v2]
                @test nin(v5) == [nout(v4), nout(v2)] == [5, 5]

                @test_logs (:warn, r"Can not add edge") create_edge!(ve, v4)
                @test inputs(v5) == [v4, v2]
                @test nin(v5) == [nout(v4), nout(v2)] == [5, 5]
            end

            @testset "PostSelectOutputs SizeInvariant" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = iv(v1, v2, name = "v4")
                v5 = av(v4, 3, name="v5")

                create_edge!(v3, v4, strategy=PostSelectOutputs(align = DecreaseBigger(), valuefun = v -> 1:nout_org(v)))

                @test inputs(v4) == [v1, v2, v3]
                @test nin(v4) == nout.([v1, v2, v3]) == [4,4,4]
                @test in_inds(op(v4)) == out_inds.(op.([v1,v2,v3])) == [1:4, 1:4, 2:5]
            end

            @testset "PostSelectOutputs SizeInvariant post align" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = iv(v1,v2, name = "v4")
                v5 = av(v4, 3, name="v5")

                create_edge!(v3, v4, strategy=PostSelectOutputs(valuefun = v -> 1:nout_org(v)))

                @test inputs(v4) == [v1, v2, v3]
                @test nin(v4) == nout.([v1, v2, v3]) == [4,4,4]
                @test in_inds(op(v4)) == out_inds.(op.([v1,v2,v3])) == [1:4, 1:4, 2:5]
            end

            @testset "PostSelectOutputs SizeInvariant add to immutable" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 5, name="v2")
                v3 = iv(v0,v1, name = "v3")
                v4 = av(v3, 3, name="v4")

                create_edge!(v2, v3, strategy=PostSelectOutputs(align= IncreaseSmaller(), valuefun = v -> 1:nout_org(v)))

                @test inputs(v3) == [v0, v1, v2]
                @test nin(v3) == nout.([v0, v1, v2]) == [3,3,3]
                @test in_inds(op(v3)) == out_inds.(op.([v1,v1,v2])) == [1:3, 1:3, 3:5]
                @test in_inds(op(v4)) == [out_inds(op(v3))] == [1:3]
            end

            @testset "PostSelectOutputs SizeStack" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = sv(v1,v2, name = "v4")
                v5 = av(v4, 3, name="v5")

                create_edge!(v3, v4, strategy=PostSelectOutputs())

                @test inputs(v4) == [v1, v2, v3]
                @test nin(v4) == nout.([v1, v2, v3]) == [3,4,5]
                @test out_inds(op(v1)) == 1:3
                @test out_inds(op(v2)) == 1:4
                @test out_inds(op(v3)) == 1:5
                @test out_inds(op(v4)) == [1:7;-ones(Int,5)]
                @test in_inds(op(v5)) == [[1:7;-ones(Int,5)]]
            end

            @testset "PostSelectOutputs fail size align" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = sv(v1,v2, name = "v4")
                v5 = av(v4, 3, name="v5")

                @test_logs (:warn, r"Could not align size") create_edge!(v3, v4, strategy=PostSelectOutputs(align = FailAlignSizeWarn()))
                @test inputs(v4) == [v1, v2]
                @test nin(v4) == nout.([v1, v2]) == [3,4]
            end

            @testset "PostSelectOutputs fail select" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = iv(v1,v2, name = "v4")
                v5 = av(v4, 3, name="v5")

                @test_logs (:warn, r"Could not align size") create_edge!(v3, v4, strategy=PostSelectOutputs(align=IncreaseSmaller(), select=NoutRevert(), fallback=FailAlignSizeWarn()))
                @test inputs(v4) == [v1, v2]
                @test nin(v4) == nout.([v1, v2]) == [4,4]
            end

            @testset "PostSelectOutputs SizeStack accidental aligned last" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = sv(v1,v2,v3, name = "v4")
                v5 = av(v4, 3, name="v5")

                # These kinda shenanigans is what one might do when one wants to remove an edge without replacing the input neurons of the next layer.
                insert!(v2, v -> conc(v, dims=1), reverse)
                dummy = outputs(v2)[]
                remove_edge!(v2, dummy; strategy = NoSizeChange())

                @test  create_edge!(v3, dummy, strategy = PostAlignJuMP())
                @test remove!(dummy, RemoveStrategy(NoSizeChange()))

                Δoutputs(v4, v -> ones(nout_org(v)))
                @test inputs(v4) == [v1, v3, v3]
                @test nin(v4) == nout.([v1, v3, v3]) == [3,5,5]
                @test nin(v5) == [nout(v4)] == [13]
            end

            @testset "PostApplyMutation SizeStack" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 3, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = sv(v1,v2, name = "v4")
                v5 = av(v4, 3, name="v5")

                create_edge!(v3, v4, strategy=PostApplyMutation())

                @test inputs(v4) == [v1, v2, v3]
                @test nin(v4) == nin_org(v4) == nout.([v1, v2, v3]) == [3,4,5]
                @test nin(v5) == nin_org(v5) == [nout(v4)] == [nout_org(v4)] == [3+4+5]
            end

            @testset "PostApplyMutation failure" begin
                v0 = inpt(3, "v0")
                v1 = av(v0, 4, name="v1")
                v2 = av(v0, 4, name="v2")
                v3 = av(v0, 5, name="v3")
                v4 = iv(v1,v2, name = "v4")
                v5 = av(v4, 3, name="v5")

                Δnout(v2, -1)
                create_edge!(v3, v4, strategy=PostApplyMutation(FailAlignSizeRevert()))

                @test inputs(v4) == [v1, v2]
                @test nin(v4) == nout.([v1, v2]) == [3,3]
                @test nin_org(v4) == nout_org.([v1, v2]) == [4,4]
            end
        end

        @testset "With size constraints" begin

            struct SizeConstraintNoDecrease
                constraint
                nodecrease
             end
            NaiveNASlib.minΔnoutfactor(c::SizeConstraintNoDecrease) = c.constraint
            NaiveNASlib.minΔninfactor(c::SizeConstraintNoDecrease) = c.constraint
            function NaiveNASlib.compconstraint!(s, c::SizeConstraintNoDecrease, data)
                model = data.model
                v = data.vertex

                fv_out = JuMP.@variable(model, integer=true)
                JuMP.@constraint(data.model, c.constraint * fv_out ==  nout(v) - data.noutdict[v])

                ins = filter(vin -> vin in keys(data.noutdict), inputs(v))

                fv_in = JuMP.@variable(model, [1:length(ins)], integer=true)
                JuMP.@constraint(data.model, [i=1:length(ins)], c.constraint * fv_in[i] == data.noutdict[ins[i]] - nin(v)[i])

                # Not directly related to size constraint. Want to always increase sizes
                if c.nodecrease
                    JuMP.@constraint(data.model, data.noutdict[data.vertex] >= nout(data.vertex))
                end
            end
            # Can't have kwarg due to https://github.com/JuliaLang/julia/issues/32350
            av(in, outsize, constr, name="avs", nodecr=true) = av(in, outsize, name=name, comp = SizeConstraintNoDecrease(constr, nodecr))

            @testset "Edge addition" begin

                @testset "Add to nout-constrained stacking" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 2, "v1")
                    v2 = av(v0, 6, 3, "v2")
                    v3 = sv(v1, name="v3")
                    v4 = av(v3, 5, 5, "v4")
                    v5 = av(v3, 7, 7, "v5")

                    @test inputs(v3) == [v1]
                    @test minΔninfactor(v3) == 70

                    create_edge!(v2, v3)
                    @test inputs(v3) == [v1, v2]

                    @test nin(v3) == [nout(v1), nout(v2)] == [28, 15]
                    @test [nout(v3)] == nin(v4) == nin(v5) == [43]
                end

                @testset "Add immutable to nout-constrained stacking" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 2, "v1")
                    v2 = sv(v1, name="v2")
                    v3 = av(v2, 5, 5, "v3")

                    @test inputs(v2) == [v1]
                    @test minΔninfactor(v2) == 10

                    create_edge!(v0, v2)
                    @test inputs(v2) == [v1, v0]

                    @test nin(v2) == [nout(v1), nout(v0)] == [10, 3]
                    @test [nout(v2)] == nin(v3) == [13]
                end

                @testset "Add nout-constrained to stacking with one immutable output" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1", false)
                    v2 = av(v0, 10, 2, "v2", false)
                    v3 = sv(v1, name="v3")
                    v4 = av(v3, 5, 5, "v4", false)
                    v5 = imu(v3, 3, name="v5")

                    @test inputs(v3) == [v1]
                    @test ismissing(minΔnoutfactor(v3))

                    create_edge!(v2, v3)
                    @test inputs(v3) == [v1, v2]

                    @test nin(v3) == [nout(v1), nout(v2)] == [2, 6]
                    @test [nout(v3)] == nin(v4) == nin(v5) == [8]
                end

                @testset "Add nout-constrained to stacking with immutable output" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1", false)
                    v2 = av(v0, 10, 2, "v2", false)
                    v3 = av(v0, 5, 5, "v3", false)
                    v4 = sv(v1,v2, name="v4")
                    v5 = imu(v4, 3, name="v5")

                    @test inputs(v4) == [v1, v2]
                    @test ismissing(minΔnoutfactor(v4))

                    create_edge!(v3, v4)
                    @test inputs(v4) == [v1, v2, v3]

                    @test nin(v4) == [nout(v1), nout(v2), nout(v3)] == [5, 8, 5]
                    @test [nout(v4)] == nin(v5) == [18]
                end

                @testset "Fail for impossible size constraint" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1")
                    v2 = av(v0, 11, 2, "v2")
                    v3 = sv(v1, name="v3")
                    v4 = imu(v3, 3, name="v4")

                    @test inputs(v3) == [v1]
                    @test ismissing(minΔnoutfactor(v3))

                    @test_throws ErrorException create_edge!(v2, v3)
                end

                @testset "Warn for impossible size constraint and revert" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1")
                    v2 = av(v0, 11, 2, "v2")
                    v3 = sv(v1, name="v3")
                    v4 = imu(v3, 3, name="v4")

                    @test inputs(v3) == [v1]
                    @test [nout(v1)] == nin(v3) == [nout(v3)] == nin(v4) == [8]
                    @test ismissing(minΔnoutfactor(v3))

                    @test_logs (:warn, r"Could not align sizes") create_edge!(v2, v3, strategy = PostAlignJuMP(DefaultJuMPΔSizeStrategy(), fallback=FailAlignSizeWarn()))

                    @test inputs(v3) == [v1]
                    @test [nout(v1)] == nin(v3) == [nout(v3)] == nin(v4) == [8]
                end

                @testset "Add to nout-constrained invariant" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 2, "v1")
                    v2 = av(v0, 6, 3, "v2")
                    v3 = iv(v1, name="v3")
                    v4 = av(v3, 5, 5, "v4")
                    v5 = av(v3, 7, 7, "v5")

                    @test inputs(v3) == [v1]
                    @test minΔninfactor(v3) == 70

                    create_edge!(v2, v3)
                    @test inputs(v3) == [v1, v2]

                    @test nin(v3) == [nout(v1), nout(v2)] == [78, 78]
                    @test [nout(v3)] == nin(v4) == nin(v5) == [78]
                end

                @testset "Add immutable to nout-constrained invariant" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 9, 3, "v1", false)
                    v2 = iv(v1, name="v2")
                    v3 = av(v2, 5, 6, "v3", false)

                    @test inputs(v2) == [v1]
                    @test minΔninfactor(v2) == 6

                    create_edge!(v0, v2)
                    @test inputs(v2) == [v1, v0]

                    @test nin(v2) == [nout(v1), nout(v0)] == [3, 3]
                    @test [nout(v2)] == nin(v3) == [3]
                end

                @testset "Add nout-constrained to invariant with one immutable output" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1", false)
                    v2 = av(v0, 13, 5, "v2", false)
                    v3 = iv(v1, name="v3")
                    v4 = av(v3, 5, 5, "v4", false)
                    v5 = imu(v3, 3, name="v5")

                    @test inputs(v3) == [v1]
                    @test ismissing(minΔnoutfactor(v3))

                    create_edge!(v2, v3)
                    @test inputs(v3) == [v1, v2]

                    @test nin(v3) == [nout(v1), nout(v2)] == [8, 8]
                    @test [nout(v3)] == nin(v4) == nin(v5) == [8]
                end

                @testset "Add nout-constrained to invariant with immutable output" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 10, 3, "v1", false)
                    v2 = av(v0, 10, 2, "v2", false)
                    v3 = av(v0, 5, 5, "v3", false)
                    v4 = iv(v1,v2, name="v4")
                    v5 = imu(v4, 3, name="v5")

                    @test inputs(v4) == [v1, v2]
                    @test ismissing(minΔnoutfactor(v4))

                    create_edge!(v3, v4)
                    @test inputs(v4) == [v1, v2, v3]

                    @test nin(v4) == [nout(v1), nout(v2), nout(v3)] == [10, 10, 10]
                    @test [nout(v4)] == nin(v5) == [10]
                end

                @testset "Fail for impossible size constraint" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1")
                    v2 = av(v0, 11, 2, "v2")
                    v3 = iv(v1, name="v3")
                    v4 = imu(v3, 3, name="v4")

                    @test inputs(v3) == [v1]
                    @test ismissing(minΔnoutfactor(v3))

                    @test_throws ErrorException create_edge!(v2, v3)
                end

                @testset "Warn for impossible size constraint and ignore" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1")
                    v2 = av(v0, 11, 2, "v2")
                    v3 = iv(v1, name="v3")
                    v4 = imu(v3, 3, name="v4")

                    @test inputs(v3) == [v1]
                    @test [nout(v1)] == nin(v3) == [nout(v3)] == nin(v4) == [8]
                    @test ismissing(minΔnoutfactor(v3))

                    @test_logs (:warn, r"Could not align sizes") create_edge!(v2, v3, strategy = PostAlignJuMP(DefaultJuMPΔSizeStrategy(), fallback=FailAlignSizeWarn()))

                    @test inputs(v3) == [v1]
                    @test [nout(v1)] == nin(v3) == [nout(v3)] == nin(v4) == [8]
                end
            end

            @testset "Edge removal" begin
                #Why no tests for Invariant here? Because constraints are moot when removing edges as everything has the same size to begin with -> no need to change anything

                @testset "Remove from nout-constrained stacking" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 5, "v1")
                    v2 = av(v0, 6, 3, "v2")
                    v3 = sv(v1, v2, name="v3")
                    v4 = av(v3, 5, 3, "v4")
                    v5 = av(v2, 7, 7, "v5")

                    @test inputs(v3) == [v1, v2]
                    @test nin(v4) == [nout(v3)] == [nout(v1) + nout(v2)] == [14]

                    remove_edge!(v2, v3)
                    @test inputs(v3) == [v1]

                    @test nin(v4) == nin(v3) == [nout(v1)] == [8]
                    @test nin(v5) == [nout(v2)] == [6] # Not touched
                end

                @testset "Remove from nout-constrained stacking with immutable input" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 2, "v1", false)
                    v2 = sv(v0, v1, name="v2")
                    v3 = av(v2, 5, 2, "v3", false)

                    @test inputs(v2) == [v0, v1]

                    remove_edge!(v1, v2)
                    @test inputs(v2) == [v0]

                    @test nin(v2) == [nout(v0)] == [3]
                    @test [nout(v2)] == nin(v3) == [3]
                end

                @testset "Remove from nout-constrained to stacking with one immutable output" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 2, "v1", false)
                    v2 = av(v0, 10, 1, "v2", false) # Δfactor of 1 is of no help as v2 will be remove
                    v3 = sv(v1,v2, name="v3")
                    v4 = av(v3, 5, 1, "v4", false) # Δfactor of 1 is of no help as v5 is immutable
                    v5 = imu(v3, 3, name="v5")

                    @test inputs(v3) == [v1,v2]
                    @test ismissing(minΔnoutfactor(v3))

                    remove_edge!(v2, v3)
                    @test inputs(v3) == [v1]

                    @test nin(v3) == [nout(v1)] == [18]
                    @test [nout(v3)] == nin(v4) == nin(v5) == [18]
                end

                @testset "Remove from nout-constrained to stacking with immutable output" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 1, "v1", false)
                    v2 = av(v0, 10, 2, "v2", false)
                    v3 = av(v0, 5, 3, "v3", false)
                    v4 = sv(v1,v2,v3, name="v4")
                    v5 = imu(v4, 3, name="v5")

                    @test inputs(v4) == [v1, v2, v3]
                    @test [nout(v4)] == nin(v5) == [23]
                    @test ismissing(minΔnoutfactor(v4))

                    remove_edge!(v1, v4)
                    @test inputs(v4) == [v2, v3]

                    @test nin(v4) == [nout(v2), nout(v3)] == [18, 5]
                    @test [nout(v4)] == nin(v5) == [23] # No change allowed
                end

                @testset "Fail for impossible size constraint" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1", false)
                    v2 = av(v0, 11, 2, "v2", false)
                    v3 = sv(v1,v2, name="v3")
                    v4 = imu(v3, 3, name="v4")

                    @test inputs(v3) == [v1, v2]
                    @test ismissing(minΔnoutfactor(v3))

                    @test_throws ErrorException remove_edge!(v2, v3)
                end

                @testset "Warn for impossible size constraint and revert" begin
                    v0 = inpt(3, "v0")
                    v1 = av(v0, 8, 3, "v1", false)
                    v2 = av(v0, 11, 2, "v2", false)
                    v3 = sv(v1, v2, name="v3")
                    v4 = imu(v3, 3, name="v4")

                    @test inputs(v3) == [v1,v2]
                    @test [nout(v1), nout(v2)] == nin(v3) == [8, 11]
                    @test [nout(v3)] == nin(v4) == [19]
                    @test ismissing(minΔnoutfactor(v3))

                    @test_logs (:warn, r"Could not align sizes") remove_edge!(v2, v3, strategy = PostAlignJuMP(DefaultJuMPΔSizeStrategy(), fallback=FailAlignSizeWarn()))

                    @test inputs(v3) == [v1,v2]
                    @test [nout(v1), nout(v2)] == nin(v3) == [8, 11]
                    @test [nout(v3)] == nin(v4) == [19]
                end
            end
        end
    end

    @testset "Vertex addition"  begin

        @testset "Add to linear graph" begin
            v0 = inpt(3, "v0")
            v1 = av(v0, 5, name="v1")
            v2 = av(v1, 4, name="v2")

            @test inputs(v2) != outputs(v1)
            graph = CompGraph(v0, v2)

            @test graph(3) == 3

            insert!(v1, v -> av(v, nout(v), name="vnew1"))

            @test inputs(v2) == outputs(v1)
            vnew1 = inputs(v2)[]
            @test [nout(v1)] == nin(vnew1) == [nout(vnew1)] == nin(v2) == [5]

            @test graph(3) == 3

            @test inputs(vnew1) == [v1]
            @test outputs(vnew1) == [v2]

            # Add two consecutive vertices
            insert!(vnew1, v -> av(av(v, 3, name="vnew2"), nout(v), name="vnew3"))

            @test [inputs(v2)] == outputs.(outputs(vnew1))
            vnew2 = outputs(vnew1)[]
            vnew3 = outputs(vnew2)[]

            @test [nout(vnew1)] == nin(vnew2) == [nout(vnew3)] == nin(v2) == [5]

            @test graph(3) == 3

        end

        @testset "Add to one of many inputs" begin
            v0 = inpt(3, "v0")
            v1 = av(v0, 4, name="v1")
            v2 = av(v0, 5, name="v2")
            v3 = av(v0, 6, name="v3")
            v4 = sv(v1,v2,v3, name="v4")
            v5 = av(v4, 7, name="v5")

            insert!(v2, v -> av(v, nout(v), name="vnew1"))
            vnew1 = outputs(v2)[]

            @test vnew1 != v2
            @test inputs(v4) == [v1, vnew1, v3]

            insert!(v0, v -> av(v, nout(v), name="vnew2"))

            @test length(outputs(v0)) == 1
            vnew2 = outputs(v0)[]

            @test inputs(v1) == inputs(v2) == inputs(v3) == [vnew2]
            @test outputs(vnew2) == [v1, v2, v3]

            insert!(vnew2, v -> av(v, nout(v), name="vnew3"), vouts -> vouts[[1, 3]])

            @test length(outputs(vnew2)) == 2
            @test outputs(vnew2)[1] == v2
            vnew3 = outputs(vnew2)[2]

            @test outputs(vnew3) == [v1, v3]
            @test inputs(v4) == [v1, vnew1, v3]
        end
    end

    @testset "Vertex removal" begin

        @testset "Fail noop" begin
            v0 = inpt(3, "v0")
            v1 = av(v0, 5, name="v1")
            v2 = av(v1, 4, name="v2")

            @test remove!(v1, RemoveStrategy(FailAlignSizeNoOp())) == false

            @test name.(inputs(v2)) == [name(v1)]
            @test name.(outputs(v1)) == [name(v2)]
        end

        @testset "Remove from linear graph" begin
            v0 = inpt(3)
            v1 = av(v0, 5, name="v1")
            v2 = av(v1, 4, name="v2")
            v3 = av(v2, 6, name="v3")

            remove!(v2)
            @test inputs(v3) == [v1]
            @test outputs(v1) == [v3]
            @test nin(v3) == [nout(v1)] == [5]

            # Note, input to v1 can not be changed, we must decrease
            # nin of v3
            remove!(v1)
            @test inputs(v3) == [v0]
            @test outputs(v0) == [v3]
            @test nin(v3) == [nout(v0)] == [3]
        end

        @testset "Remove one of many inputs to stacking" begin
            v0 = inpt(3, "v0")
            v1 = av(v0, 4, name="v1")
            v2 = av(v0, 5, name="v2")
            v3 = av(v0, 6, name="v3")
            v4 = sv(v1,v2,v3, name="v4")
            v5 = av(v4, 7, name="v5")

            # Note, input to v1 can not be changed, we must decrease
            # nin of v4 (and v5)
            remove!(v2)
            @test inputs(v4) == [v1, v0, v3]
            @test nin(v4) == [nout(v1), nout(v0), nout(v3)] == [4,3,6]
            @test nin(v5) == [nout(v4)] == [3+4+6]

            #Now lets try without connecting the inputs to v4
            remove!(v1, RemoveStrategy(ConnectNone(), ChangeNinOfOutputs(-nout(v1))))
            @test inputs(v4) == [v0, v3]
            @test nin(v4) == [nout(v0), nout(v3)] == [3, 6]
            @test nin(v5) == [nout(v4)] == [3+6]
        end

        @testset "Remove one of many inputs to stacking increase size" begin
            v0 = inpt(7, "v0")
            v1 = av(v0, 4, name="v1")
            v2 = av(v0, 5, name="v2")
            v3 = av(v0, 6, name="v3")
            v4 = sv(v1,v2,v3, name="v4")
            v5 = av(v4, 7, name="v5")

            # Note, input to v1 can not be changed, we must decrease
            # nin of v4 (and v5)
            remove!(v2)
            @test inputs(v4) == [v1, v0, v3]
            @test nin(v4) == [nout(v1), nout(v0), nout(v3)] == [4,7,6]
            @test nin(v5) == [nout(v4)] == [4+7+6]

            #Now lets try without connecting the inputs to v4
            remove!(v1, RemoveStrategy(ConnectNone(), ChangeNinOfOutputs(-nout(v1))))
            @test inputs(v4) == [v0, v3]
            @test nin(v4) == [nout(v0), nout(v3)] == [7, 6]
            @test nin(v5) == [nout(v4)] == [7+6]
        end

        @testset "Remove one of many inputs to invariant" begin
            v0 = inpt(3, "v0")
            v1 = av(v0, 5, name="v1")
            v2 = av(v0, 5, name="v2")
            v3 = av(v0, 5, name="v3")
            v4 = iv(v1,v2,v3, name="v4")
            v5 = av(v4, 7, name="v5")

            # Note, input to v1 can not be changed, we must decrease
            # nin of v4 (and v5)
            remove!(v2)
            @test inputs(v4) == [v1, v0, v3]
            @test nin(v4) == [nout(v1), nout(v0), nout(v3)] == [3,3,3]
            @test nin(v5) == [nout(v4)] == [3]

            #Now lets try without connecting the inputs to v4
            # NoSizeChange is just to avoid touching the input vertex
            remove!(v1, RemoveStrategy(ConnectNone(), NoSizeChange()))
            @test inputs(v4) == [v0, v3]
            @test nin(v4) == [nout(v0), nout(v3)] == [3, 3]
            @test nin(v5) == [nout(v4)] == [3]
        end

        @testset "Remove one of many inputs to invariant increase size" begin
            v0 = inpt(7, "v0")
            v1 = av(v0, 5, name="v1")
            v2 = av(v0, 5, name="v2")
            v3 = av(v0, 5, name="v3")
            v4 = iv(v1,v2,v3, name="v4")
            v5 = av(v4, 7, name="v5")

            # Note, input to v1 can not be changed, we must decrease
            # nin of v4 (and v5)
            remove!(v2)
            @test inputs(v4) == [v1, v0, v3]
            @test nin(v4) == [nout(v1), nout(v0), nout(v3)] == [7,7,7]
            @test nin(v5) == [nout(v4)] == [7]

            #Now lets try without connecting the inputs to v4
            # NoSizeChange is just to avoid touching the input vertex
            remove!(v1, RemoveStrategy(ConnectNone(), NoSizeChange()))
            @test inputs(v4) == [v0, v3]
            @test nin(v4) == [nout(v0), nout(v3)] == [7, 7]
            @test nin(v5) == [nout(v4)] == [7]
        end

        @testset "Remove input duplicated stacking" begin
            v0 = inpt(3, "v0")
            v1 = av(v0, 4, name="v1")
            v2 = av(v0, 5, name="v2")
            v3 = sv(v1,v2,v2,v1, name="v3")
            v4 = av(v3, 7, name="v4")

            remove!(v1)
            @test inputs(v3) == [v0, v2, v2, v0]
            @test nin(v3) == [nout(v0), nout(v2), nout(v2), nout(v0)] == [3,5,5,3]
            @test nin(v4) == [nout(v3)] == [3+5+5+3]

            #Now lets try without connecting the inputs to v3
            remove!(v2, RemoveStrategy(ConnectNone(), ChangeNinOfOutputs(-nout(v2))))
            @test inputs(v3) == [v0, v0]
            @test nin(v3) == [nout(v0), nout(v0)] == [3,3]
            @test nin(v4) == [nout(v3)] == [3+3]
        end

        @testset "Remove input duplicated invariant" begin
            v0 = inpt(3, "v0")
            v1 = av(v0, 5, name="v1")
            v2 = av(v0, 5, name="v2")
            v3 = iv(v1,v2,v2,v1, name="v3")
            v4 = av(v3, 7, name="v4")

            remove!(v1)
            @test inputs(v3) == [v0, v2, v2, v0]
            @test nin(v3) == [nout(v0), nout(v2), nout(v2), nout(v0)] == [3,3,3,3]
            @test nin(v4) == [nout(v3)] == [3]

            #Now lets try without connecting the inputs to v3
            # NoSizeChange is just to avoid touching the input vertex
            remove!(v2, RemoveStrategy(ConnectNone(), NoSizeChange()))
            @test inputs(v3) == [v0, v0]
            @test nin(v3) == [nout(v0), nout(v0)] == [3,3]
            @test nin(v4) == [nout(v3)] == [3]
        end

        @testset "Remove one of many outputs" begin
            v0 = inpt(3)
            v1 = av(v0, 4)
            v2 = av(v1, 5)
            v3 = av(v1, 6)
            v4 = av(v1, 7)
            v5 = av(v2, 8)

            remove!(v2)
            @test outputs(v1) == [v5, v3, v4]
            @test nin(v5) == nin(v3) == nin(v4) == [nout(v1)] == [5]

            # Test that it is possible to remove vertex without any outputs
            remove!(v3)
            @test outputs(v1) == [v5, v4]
            @test nin(v5) == nin(v4) == [nout(v1)] == [6]
        end

        @testset "Remove one of many outputs PostAlignJuMP" begin
            v0 = inpt(3)
            v1 = av(v0, 4, name="v1")
            v2 = av(v1, 5, name="v2")
            v3 = av(v1, 6, name="v3")
            v4 = av(v1, 7, name="v4")
            v5 = av(v2, 8, name="v5")

            remove!(v2, RemoveStrategy(PostAlignJuMP()))
            @test outputs(v1) == [v5, v3, v4]
            @test nin(v5) == nin(v3) == nin(v4) == [nout(v1)] == [4]

            # Test that it is possible to remove vertex without any outputs
            remove!(v3, RemoveStrategy(PostAlignJuMP()))
            @test outputs(v1) == [v5, v4]
            @test nin(v5) == nin(v4) == [nout(v1)] == [4]
        end

        @testset "Hidden immutable" begin
            v0 = inpt(3)
            v1 = sv(v0)
            v2 = av(v1, 4)
            v3 = av(v2, 5)

            #Danger! Must realize that size of v1 can not be changed!
            remove!(v2)
            @test outputs(v1) == [v3]
            @test inputs(v3) == [v1]
            @test nin(v3) == [nout(v1)] == [3]
        end

        @testset "Incompatible size factor" begin
            v1 = av(inpt(3), 5, name="v1")
            p1 = iv(v1, name="p1")
            p2 = iv(v1, name="p2")
            p3 = iv(v1, name="p3")
            join = sv(p1,p2,p3, name="join")
            v2 = av(join, 16, name = "v2") # 16 is not divisible by 3!
            v3 = av(v2, 4, name="v3")

            @test minΔnoutfactor_only_for.(outputs(v2)) == [1]
            @test minΔnoutfactor_only_for.(inputs(v2)) == [3]

            # Impossible to set nout of join to 16 as it is a join of the same vertex 3 times (which is obviously a senseless construct)
            remove!(v2)
            @test nin(v3) == [nout(join)] == [3nout(v1)] == [15]
        end

        @testset "Fail for impossible removal" begin
            v1 = inpt(3, "v1")
            v2 = av(v1, 5, name="v2")
            v3 = imu(v2, 4, name="v3")

            @test_throws ErrorException remove!(v2)
        end

        @testset "Warn for impossible removal and ignore" begin
            v1 = inpt(3, "v1")
            v2 = av(v1, 5, name="v2")
            v3 = imu(v2, 4, name="v3")

            @test_logs (:warn, r"Could not align sizes") remove!(v2, RemoveStrategy(AlignSizeBoth(FailAlignSizeWarn())))

            @test outputs(v1) == [v2]
            @test inputs(v2) == [v1]
            @test outputs(v2) == [v3]
            @test inputs(v3) == [v2]
            @test [nout(v1)] == nin(v2) == [3]
            @test [nout(v2)] == nin(v3) == [5]
        end

        @testset "Size constraints" begin

            struct SizeConstraint constraint; end
            NaiveNASlib.minΔnoutfactor(c::SizeConstraint) = c.constraint
            NaiveNASlib.minΔninfactor(c::SizeConstraint) = c.constraint

            @testset "Incompatible size constraints" begin

                v1 = av(inpt(3), 10, name="v1", comp = SizeConstraint(2))
                v2 = av(v1, 5, name = "v2")
                v3 = av(v2, 4, name="v3", comp = SizeConstraint(3))

                @test minΔnoutfactor_only_for.(outputs(v2)) == [3]
                @test minΔnoutfactor_only_for.(inputs(v2)) == [2]

                # Impossible to increase v1 by 5 due to SizeConstraint(3)
                # But also impossible to decrease nin of v3 by 5 due to SizeConstraint(2)
                # However, if we decrease v1 by 2 and increase v3 by 3 we will hit home!
                # Fallback to AlignBoth which does just that
                remove!(v2)
                @test nin(v3) == [nout(v1)] == [8]
            end

            @testset "Incompatible size constraints transparent vertex" begin

                v1 = av(inpt(3), 10, name="v1", comp = SizeConstraint(2))
                v2 = sv(v1, name = "v2")
                v3 = av(v2, 4, name="v3", comp = SizeConstraint(3))

                @test minΔnoutfactor_only_for.(outputs(v2)) == [3]
                @test minΔnoutfactor_only_for.(inputs(v2)) == [2]

                # Size is already aligned due to transparent. Just test that this
                # does not muck things up
                remove!(v2, RemoveStrategy(AlignSizeBoth()))
                @test nin(v3) == [nout(v1)] == [10]
            end
        end

        @testset "Tricky structures" begin

            @testset "Remove residual layers" begin
                v1 = av(inpt(3, "in"), 10, name="v1")
                v2 = av(v1, 3, name="v2")
                v3 = sv(v2, name="v3")
                v4 = av(v3, 10, name="v4")
                v5 = iv(v4, v1, name="v5")
                v6 = av(v5, 4, name="v6")

                remove!(v4)
                @test inputs(v5) == [v3, v1]
                @test nin(v5) == [nout(v3), nout(v1)] == [10, 10]
                @test nin(v6) == [nout(v5)] == [10]

                remove!(v3)
                @test inputs(v5) == [v2, v1]
                @test nin(v5) == [nout(v2), nout(v1)] == [10, 10]
                @test nin(v6) == [nout(v5)] == [10]

                remove!(v2)
                @test inputs(v5) == [v1, v1]
                @test nin(v5) == [nout(v1), nout(v1)] == [10, 10]
                @test nin(v6) == [nout(v5)] == [10]

                v7 = av(v6, 13, name="v7")
                remove!(v6)
                @test nin(v5) == [nout(v1), nout(v1)] == [10, 10]
                @test nin(v7) == [nout(v5)] == [10]

                v8 = av(v7, 3, name="v8")
                remove!(v7)
                @test nin(v5) == [nout(v1), nout(v1)] == [13, 13]
                @test nin(v7) == [nout(v5)] == [13]
            end

            @testset "Remove after transparent fork" begin
                v1 = av(inpt(3, "in"), 5, name="v1")
                p1 = iv(v1, name="p1")
                p2 = iv(v1, name="p2")
                v2 = sv(p1,p2, name="v2")
                v3 = iv(v2, name="v3")
                v4 = av(v3, 12, name="v4")
                v5 = av(v4, 7, name="v5")

                remove!(v4)
                @test inputs(v5) == [v3]
                @test outputs(v3) == [v5]

                @test nin(v5) == [nout(v3)] == [nout(p1) + nout(p2)] == [2 * nout(v1)] == [12]
            end

            @testset "Remove before half transparent resblock" begin
                v1 = av(inpt(2, "in"), 5, name="v1")
                v2 = av(v1, 3, name="v2")
                v3 = iv(v2, name="v3")
                v4 = av(v3, 3, name="v4")
                v5 = iv(v4, name="v5")
                v6 = iv(v3, v4, name="v6")
                v7 = av(v6, 2, name="v7")

                remove!(v5)
                @test inputs(v6) == [v3, v4]
                @test nin(v6) == [nout(v3), nout(v4)] == [3, 3]
            end

            @testset "Remove right before fork" begin
                v1 = av(inpt(3, "in"), 3, name="v1")
                v2 = av(v1, 5, name="v2")
                p1 = iv(v2, name="p1")
                p2₁ = iv(v2, name="p2_1")
                p2₂ = av(p2₁,7, name="p2_2")
                p2₃ = iv(p2₂, name="p2_3")
                v3 = sv(p1, p2₃, name="v3")
                v4 = iv(v3, name="v4")
                v5 = av(v4, 12, name="v5")

                remove!(v2)
                @test inputs(p1) == inputs(p2₁) == [v1]
                @test nin(p1) == nin(p2₁) == [nout(v1)] == [5]
            end

            @testset "Remove nout cycle back to nin aligned" begin
                v1 = av(inpt(3, "in"), 3, name="v1")
                v2 = av(v1, 5, name="v2")
                v3 = iv(v2, name="v3")
                v4 = av(v3, 5, name="v4")
                v5 = iv(v4, name="v5")
                v6 = "v6" >> v5 + v2

                remove!(v4, RemoveStrategy(CheckNoSizeCycle()))
                @test inputs(v5) == [v3]
                @test nin(v5) == [nout(v3)] == [5]
            end

            @testset "Remove nout cycle back to nin abort" begin
                v1 = av(inpt(3, "in"), 3, name="v1")
                v2 = av(v1, 5, name="v2")
                v3 = iv(v2, name="v3")
                p1 = av(v3, 3, name="p1")
                p2 = av(v3, 2, name="p2")
                v4 = sv(p1,p2, name="v4")
                v5 = "v5" >> v4 + v2

                # Evilness: Removing p1 would create a "size transparent loop" where nout(v4) = nout(v4) + nout(p2)
                @test_logs (:warn, "Can not remove vertex $(p1)! Size cycle detected!") remove!(p1)
                @test inputs(v4) == [p1, p2]
                @test nout.(inputs(v4)) == nin(v4) == [3,2]
                @test sum(nin(v4)) == nout(v4) == 5
            end

            @testset "Remove nout almost cycle back to nin" begin
                v1 = av(inpt(3, "in"), 3, name="v1")
                v2 = av(v1, 5, name="v2")
                v3 = iv(v2, name="v3")
                p1a = av(v3, 2, name="p1a")
                p1b = av(p1a, 3, name="p1b")
                p2 = av(v3, 2, name="p2")
                v4 = sv(p1b,p2, name="v4")
                v5 = "v5" >> v4 + v2

                remove!(p1b)
                @test inputs(v4) == [p1a, p2]
                @test nout.(inputs(v4)) == nin(v4) == [3,2]
                @test sum(nin(v4)) == nout(v4) == 5
            end

            @testset "Remove nout almost cycle back to nin" begin
                v1 = av(inpt(3, "in"), 3, name="v1")
                v2 = av(v1, 5, name="v2")
                v3 = iv(v2, name="v3")
                p1a = av(v3, 2, name="p1a")
                p1b = av(p1a, 3, name="p1b")
                p2 = av(v3, 2, name="p2")
                v4 = sv(p1b,p2, name="v4")
                v5 = "v5" >> v4 + v2

                remove!(p1a)
                @test inputs(v4) == [p1b, p2]
                @test nout.(inputs(v4)) == nin(v4) == [3,2]
                @test sum(nin(v4)) == nout(v4) == 5
            end
        end
    end
end

ext: jl


text:
module PushballNN

import SearchBall, GeometryTypes, Chipmunk, LinearAlgebra
const sb = SearchBall
const gt = GeometryTypes
const Cp = Chipmunk
const la = LinearAlgebra

# function init_strategies(config, w=rand(2,5).*2.-1, b=rand(2).*2.-1)
function init_strategies(config; kwargs...)
    x = haskey(kwargs, :x) ? kwargs[:x] : rand(12).*2 .- 1
    issimulation = haskey(kwargs, :issimulation) ? kwargs[:issimulation] : false

    w = reshape(@view(x[1:10]), (2, 5))
    b = @view x[11:12]
    # right goal
    rightgoalbody = Cp.StaticBody()
    rightgoalbackshape = Cp.SegmentShape(rightgoalbody, Cp.Vect(9.7, 2.25), Cp.Vect(9.7, 3.75), 0.15)
    Cp.set_friction(rightgoalbackshape, 1)
    # Cp.set_collision_type(rightgoalbackshape, 2)
    Cp.add_shape(config.space, rightgoalbackshape)
    upperrightgoalshape = Cp.SegmentShape(rightgoalbody, Cp.Vect(9.0, 2.25), Cp.Vect(9.7, 2.25), 0.05)
    Cp.set_friction(upperrightgoalshape, 1)
    # Cp.set_collision_type(upperrightgoalshape, 2)
    Cp.add_shape(config.space, upperrightgoalshape)
    lowerrightgoalshape = Cp.SegmentShape(rightgoalbody, Cp.Vect(9.0, 3.75), Cp.Vect(9.7, 3.75), 0.05)
    Cp.set_friction(lowerrightgoalshape, 1)
    # Cp.set_collision_type(lowerrightgoalshape, 2)
    Cp.add_shape(config.space, lowerrightgoalshape)
    Cp.add_body(config.space, rightgoalbody)
    # rightgoalline
    rightgoallinebody = Cp.StaticBody()
    rightgoallineshape = Cp.SegmentShape(rightgoallinebody, Cp.Vect(9.1, 2.25), Cp.Vect(9.1, 3.75), 0.0)
    Cp.set_sensor(rightgoallineshape, true)
    Cp.set_collision_type(rightgoallineshape, 3)
    Cp.add_shape(config.space, rightgoallineshape)
    Cp.add_body(config.space, rightgoallinebody)

    # left goal
    leftgoal_body = Cp.StaticBody()
    leftgoalbackshape = Cp.SegmentShape(leftgoal_body, Cp.Vect(-0.7, 2.25), Cp.Vect(-0.7, 3.75), 0.15)
    Cp.set_friction(leftgoalbackshape, 1)
    # Cp.set_collision_type(leftgoalbackshape, 2)
    Cp.add_shape(config.space, leftgoalbackshape)
    upperleftgoalshape = Cp.SegmentShape(leftgoal_body, Cp.Vect(0.0, 2.25), Cp.Vect(-0.7, 2.25), 0.05)
    Cp.set_friction(upperleftgoalshape, 1)
    # Cp.set_collision_type(upperleftgoalshape, 2)
    Cp.add_shape(config.space, upperleftgoalshape)
    lowerleftgoalshape = Cp.SegmentShape(leftgoal_body, Cp.Vect(0.0, 3.75), Cp.Vect(-0.7, 3.75), 0.05)
    Cp.set_friction(lowerleftgoalshape, 1)
    # Cp.set_collision_type(lowerleftgoalshape, 2)
    Cp.add_shape(config.space, lowerleftgoalshape)
    Cp.add_body(config.space, leftgoal_body)
    # leftgoalline
    leftgoallinebody = Cp.StaticBody()
    leftgoallineshape = Cp.SegmentShape(leftgoallinebody, Cp.Vect(-0.1, 2.25), Cp.Vect(-0.1, 3.75), 0.0)
    Cp.set_sensor(leftgoallineshape, true)
    Cp.set_collision_type(leftgoallineshape, 3)
    Cp.add_shape(config.space, leftgoallineshape)
    Cp.add_body(config.space, leftgoallinebody)

    if !issimulation
        sb.show_graphic!(config, "rightgoal1", [[9.7, 2.25], [9.7, 3.75]], sb.GraphicLine(linewidth=10.0, color=[0.0, 0.0, 0.0, 1.0]))
        sb.show_graphic!(config, "rightgoal2", [[9.0, 2.25], [9.7, 2.25]], sb.GraphicLine(linewidth=10.0, color=[0.0, 0.0, 0.0, 1.0]))
        sb.show_graphic!(config, "rightgoal3", [[9.0, 3.75], [9.7, 3.75]], sb.GraphicLine(linewidth=10.0, color=[0.0, 0.0, 0.0, 1.0]))
        sb.show_graphic!(config, "leftgoal1", [[-0.7, 2.25], [-0.7, 3.75]], sb.GraphicLine(linewidth=10.0, color=[0.0, 0.0, 0.0, 1.0]))
        sb.show_graphic!(config, "leftgoal2", [[0.0, 2.25], [-0.7, 2.25]], sb.GraphicLine(linewidth=10.0, color=[0.0, 0.0, 0.0, 1.0]))
        sb.show_graphic!(config, "leftgoal3", [[0.0, 3.75], [-0.7, 3.75]], sb.GraphicLine(linewidth=10.0, color=[0.0, 0.0, 0.0, 1.0]))

        sb.init_text!(config, "score", 0.0)
        sb.init_text!(config, "goal", false)
        sb.init_text!(config, "antigoal", false)
        sb.init_text!(config, "out", false)

        # button for show vectorfield
        sb.add_check_button!(config, "show Vectorfield", "show Vectorfield") do widget
            run_strategies(config)
        end
    end

    outlinebody = Cp.StaticBody()
    upoutlineshape = Cp.SegmentShape(outlinebody, Cp.Vect(-0.1, 6.1), Cp.Vect(9.1, 6.1), 0.0)
    Cp.set_sensor(upoutlineshape, true)
    Cp.set_collision_type(upoutlineshape, 4)
    Cp.add_shape(config.space, upoutlineshape)
    lowoutlineshape = Cp.SegmentShape(outlinebody, Cp.Vect(-0.1, -0.1), Cp.Vect(9.1, -0.1), 0.0)
    Cp.set_sensor(lowoutlineshape, true)
    Cp.set_collision_type(lowoutlineshape, 4)
    Cp.add_shape(config.space, lowoutlineshape)
    Cp.add_body(config.space, outlinebody)

    goaloutlinebody = Cp.StaticBody()
    upgoaloutlineshape_left = Cp.SegmentShape(goaloutlinebody, Cp.Vect(-0.1, 6.1), Cp.Vect(-0.1, 3.75), 0.0)
    Cp.set_sensor(upgoaloutlineshape_left, true)
    Cp.set_collision_type(upgoaloutlineshape_left, 5)
    Cp.add_shape(config.space, upgoaloutlineshape_left)
    lowgoaloutlineshape_left = Cp.SegmentShape(goaloutlinebody, Cp.Vect(-0.1, 2.25), Cp.Vect(-0.1, -0.1), 0.0)
    Cp.set_sensor(lowgoaloutlineshape_left, true)
    Cp.set_collision_type(lowgoaloutlineshape_left, 5)
    Cp.add_shape(config.space, lowgoaloutlineshape_left)
    upgoaloutlineshape_right = Cp.SegmentShape(goaloutlinebody, Cp.Vect(9.1, 6.1), Cp.Vect(9.1, 3.75), 0.0)
    Cp.set_sensor(upgoaloutlineshape_right, true)
    Cp.set_collision_type(upgoaloutlineshape_right, 5)
    Cp.add_shape(config.space, upgoaloutlineshape_right)
    lowgoaloutlineshape_right = Cp.SegmentShape(goaloutlinebody, Cp.Vect(9.1, 2.25), Cp.Vect(9.1, -0.1), 0.0)
    Cp.set_sensor(lowgoaloutlineshape_right, true)
    Cp.set_collision_type(lowgoaloutlineshape_right, 5)
    Cp.add_shape(config.space, lowgoaloutlineshape_right)
    Cp.add_body(config.space, goaloutlinebody)

    function mycollisionhandler(arbiter_ptr::Ptr{Nothing}, space_ptr::Ptr{Nothing}, data_ptr::Ptr{Nothing})
        arbiter = Cp.Arbiter(arbiter_ptr)
        my_config = unsafe_pointer_to_objref(data_ptr)
        if my_config.ws.obstacles[1][1] > 9.0
            my_config.global_communication["goal"] = true
            if !my_config.global_communication["issimulation"]
                sb.show_text!(my_config, "goal", true)
            end
        elseif my_config.ws.obstacles[1][1] < 0.0
            my_config.global_communication["antigoal"] = true
            if !my_config.global_communication["issimulation"]
                sb.show_text!(my_config, "antigoal", true)
            end
        end
        return true
    end
    function ballout(arbiter_ptr::Ptr{Nothing}, space_ptr::Ptr{Nothing}, data_ptr::Ptr{Nothing})
        arbiter = Cp.Arbiter(arbiter_ptr)
        space = Cp.Space(space_ptr)
        my_config = unsafe_pointer_to_objref(data_ptr)
        my_config.global_communication["out"] = true
        if !my_config.global_communication["issimulation"]
            sb.show_text!(my_config, "out", true)
        end
        return true
    end

    Cp.cpSetCollisionHandler(config.space, 1, 3, mycollisionhandler, nothing, nothing, nothing, config)
    Cp.cpSetCollisionHandler(config.space, 1, 4, ballout, nothing, nothing, nothing, config)
    Cp.cpSetCollisionHandler(config.space, 1, 5, ballout, nothing, nothing, nothing, config)

    # @show w,b
    config.global_communication["neuralnetwork"] = (dat, pos, x, y, z, m) ->predict!(dat, pos, x, y, z, m, w, b)
    config.global_communication["goal"] = false
    config.global_communication["antigoal"] = false
    config.global_communication["out"] = false
    config.global_communication["issimulation"] = issimulation
    config.global_communication["vel"] = Vector{Float64}(undef, 2)
    config.global_communication["pos"] = Vector{Float64}(undef, 5)
    config.global_communication["matprod"] = Vector{Float64}(undef, 2)

    return
end

const goal_line = [9.1, 3.0]
function objective(config::sb.Config)
    ws = config.ws
    obs = ws.obstacles[1]
    player_ball = sb.distance(ws.myPlayers[1], obs)
    ball_goal = sb.distance(obs, goal_line)
    return player_ball + 3*ball_goal
end

function predict!(res::AbstractVector{<:Number}, pos::AbstractVector{<:Number}, pos_player::AbstractVector{<:Number}, pos_ball::AbstractVector{<:Number}, x_goal::Number, matprod::AbstractVector{<:Number}, w::AbstractMatrix{<:Number}, b::AbstractVector{<:Number})
    # pos = Vector{typeof(x_goal)}(undef, 5)
    @inbounds pos[1] = (pos_player[1] - 4.5)/9.0
    @inbounds pos[2] = (pos_player[2] - 3.0)/6.0
    @inbounds pos[3] = (pos_ball[1] - 4.5)/9.0
    @inbounds pos[4] = (pos_ball[2] - 3.0)/6.0
    @inbounds pos[5] = (x_goal - 4.5)/9.0
    la.mul!(matprod, w, pos)
    res .= tanh.(matprod .+ b).*0.22 # TODO: tooks too much memory
    return
end

function calc_myPlayers(config::sb.Config, dt::Number, set_type::Symbol, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    ws = config.ws

    pos_p = ws.myPlayers[1]
    pos_b = ws.obstacles[1]

    neuralnet! = config.global_communication["neuralnetwork"]::Function
    x_g = 9.0

    vel = config.global_communication["vel"]::Vector{Float64}
    pos = config.global_communication["pos"]::Vector{Float64}
    matprod = config.global_communication["matprod"]::Vector{Float64}
    neuralnet!(vel, pos, pos_p, pos_b, x_g, matprod)

    if !config.global_communication["issimulation"]
        score = objective(config)
        sb.show_text!(config, "score", score)
        if sb.get_manip_value(config, "show Vectorfield")
            for x in sb.raster(-config.field_size[1]/2, config.field_size[1]/2, 18), y in sb.raster(-config.field_size[2]/2, config.field_size[2]/2, 12)
                veltemp = Vector{Float64}(undef, 2)
                postemp = [x, y]
                neuralnet!(veltemp, pos, postemp, pos_b, x_g)
                sb.show_graphic!(config, "arrow"*string(x, y), [postemp, postemp.+veltemp], sb.GraphicArrow(arrowloc=1.0, color=[0.0, 0.8, 0.2, 0.6], linewidth=3.0, arrowsize=0.35))
            end
        else
            sb.remove_graphics_filter!(config) do s
                startswith(s, "arrow")
            end
        end
    end

    sb.send_order!(ws.myPlayers[1], sb.SetVelocity(vel[1], vel[2])) # TODO: tooks too many memory
    # ws.myPlayers[1]
    # vell = sqrt(vel[1]^2 + vel[2]^2)
    # vell = pythag(vel)
    # z = zero(ws.myPlayers[1].max_velocity)
    # v = clamp(vell, 0.0, ws.myPlayers[1].max_velocity)/vell
    # v = clamp(vell, 0.0, 0.3)/vell
    # sb.set_velocity!(ws.myPlayers[1], v*vel[1], v*vel[2]) # TODO: tooks too many memory

    return
end

function pythag(arr::AbstractVector{T})::T where T<:Number
    return sqrt(arr[1]^2 + arr[2]^2)
end
function calc_opponend(config::sb.Config, dt::Number, set_type::Symbol, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    ws = config.ws

    # do stuff

    return
end

function run_strategies(config::sb.Config, dt::Number=0, set_type::Symbol=:default, moved_object::Union{sb.Object{<:Number}, Nothing}=nothing)
    calc_myPlayers(config, dt, set_type, moved_object)
    calc_opponend(config, dt, set_type, moved_object)
    return
end

end


text:
type SparkConf
    jconf::JSparkConf
end

function SparkConf(;opts...)
    jconf = JSparkConf(())
    opts = Dict(opts)
    for (k, v) in opts
        jcall(jconf, "set", JSparkConf, (JString, JString), string(k), v)
    end
    return SparkConf(jconf)
end

function Base.show(io::IO, conf::SparkConf)
    print(io, "SparkConf()")
end

function Base.setindex!(conf::SparkConf, key::AbstractString, val::AbstractString)
    jcall(conf.jconf, "set", JSparkConf, (JString, JString), key, val)
end

function Base.getindex(conf::SparkConf, key::AbstractString)
    jcall(conf.jconf, "get", JString, (JString,), key)
end

function Base.get(conf::SparkConf, key::AbstractString, default::AbstractString)
    jcall(conf.jconf, "get", JString, (JString, JString), key, default)
end

function setmaster(conf::SparkConf, master::AbstractString)
    jcall(conf.jconf, "setMaster", JSparkConf, (JString,), master)
end

function setappname(conf::SparkConf, appname::AbstractString)
    jcall(conf.jconf, "setAppName", JSparkConf, (JString,), appname)
end


text:
export @linq, linq

##############################################################################
##
## @linq - general macro that creates a mini DSL for chaining and macro and
##         function calls
##
##############################################################################

macro linq(arg)
    esc(replacefuns(replacechains(arg)))
end

# Snippet from Calculus.jl
type SymbolParameter{T} end
SymbolParameter(s::Symbol) = SymbolParameter{s}()

replacefuns(x) = x  # default for non-expression stuff
function replacefuns(e::Expr)
    for i in 1:length(e.args)
        e.args[i] = replacefuns(e.args[i])
    end
    if e.head == :call
        return linq(SymbolParameter(e.args[1]), e.args[2:end]...)
    else
        return e
    end
end

replacechains(x) = x
function replacechains(e::Expr)
    for i in 1:length(e.args)
        e.args[i] = replacechains(e.args[i])
    end
    if e.head == :call && e.args[1] == :|> && isa(e.args[3], Expr)
        newe = e.args[3]
        insert!(newe.args, 2, e.args[2])
        return newe
    else
        return e
    end
end

##############################################################################
##
## Various linq helper definitions
##
##############################################################################

## Default, no-op:
linq{s}(::SymbolParameter{s}, args...) = Expr(:call, s, args...)

function linq(::SymbolParameter{:with}, d, body)
    with_helper(d, body)
end

function linq(::SymbolParameter{:ix}, d, args...)
    ix_helper(d, args...)
end

function linq(::SymbolParameter{:where}, d, arg)
    where_helper(d, arg)
end

function linq(::SymbolParameter{:orderby}, d, args...)
    :(let _D = $d;  orderby(_D, _DF -> @with(_DF, DataFramesMeta.orderbyconstructor(_D)($(args...)))); end)
end

function linq(::SymbolParameter{:transform}, x, args...)
    transform_helper(x, args...)
end

function linq(::SymbolParameter{:based_on}, x, args...)
    :( DataFrames.based_on($x, _DF -> @with(_DF, DataFrame($(args...)))) )
end

function linq(::SymbolParameter{:by}, x, what, args...)
    :( by($x, $what, _DF -> @with(_DF, DataFrame($(args...)))) )
end

function linq(::SymbolParameter{:select}, x, args...)
    :(let _DF = $x; @with(_DF, select(_DF, $(expandargs(args)...))); end)
end


text:
import LinearAlgebra: ldiv!, \

export Identity

#### Type-handling
"""
    Adivtype(A, b)
Determine type of the division of an element of `b` against an element of `A`:
`typeof(one(eltype(b))/one(eltype(A)))`
"""
Adivtype(A, b) = typeof(one(eltype(b))/one(eltype(A)))

"""
    zerox(A, b)

Build a zeros vector similar to `b` of eltype `Adivtype(A,b)`.
"""
function zerox(A, b)
    T = Adivtype(A, b)
    x = similar(b, T, size(A, 2))
    fill!(x, zero(T))
    return x
end

"""
No-op preconditioner
"""
struct Identity end

\(::Identity, x) = copy(x)
ldiv!(::Identity, x) = x
ldiv!(y, ::Identity, x) = copyto!(y, x)


text:
# Export datatypes
export ELFDynEntry, ELFDynEntry32, ELFDynEntry64, ELFDynEntryRef, ELFDynEntries

# Export methods
export dyn_entry_type, dyn_entry_is_string, dyn_entry_type_string

import Base: show

"""
    ELFDynEntry

ELF Dynamic table entry type.  This is comprised by a `d_tag` member and a
`d_val` member; the `d_tag` tells what kind of command this is to the dynamic
linker (e.g. `DT_NEEDED` denotes a shared library that must be loaded for this
ELF object to link properly), whereas `d_val` is a pointer to another data
structure that contains more information.  (E.g. for a `DT_NEEDED` entry,
`d_val` would represent an offset within the dynamic string table)
"""
abstract type ELFDynEntry{H <: ELFHandle} end

@io struct ELFDynEntry32{H <: ELFHandle} <: ELFDynEntry{H}
    d_tag::UInt32
    # We drop the `d_un` union, and use only `d_val`, omitting `d_ptr`.
    d_val::UInt32
end

@io struct ELFDynEntry64{H <: ELFHandle} <: ELFDynEntry{H}
    d_tag::UInt64
    # We drop the `d_un` union, and use only `d_val`, omitting `d_ptr`.
    d_val::UInt64
end

"""
    dyn_entry_type(d::ELFDynEntry)

Return the type of the given `ELFDynEntry`
"""
dyn_entry_type(d::ELFDynEntry) = d.d_tag

"""
    dyn_entry_is_string(d::ELFDynEntry)

Return `true` if the given `ELFDynEntry` represents an offset within the
dynamic string table, and therefore can be used in a `strtab_lookup()`
"""
function dyn_entry_is_string(d::ELFDynEntry)
    string_types = [
        DT_NEEDED,
        DT_SONAME,
        DT_RPATH,
        DT_RUNPATH,
        DT_AUXILIARY,
        DT_FILTER,
        DT_CONFIG,
        DT_DEPAUDIT,
        DT_AUDIT
    ]
    return dyn_entry_type(d) in string_types
end

"""
    dyn_entry_type_string(d::ELFDynEntry)

Return the given `ELFDynEntry`'s type as a string.
"""
function dyn_entry_type_string(d::ELFDynEntry)
    global DYNAMIC_TYPE

    d_type = dyn_entry_type(d)
    if haskey(DYNAMIC_TYPE, d_type)
        return DYNAMIC_TYPE[d_type]
    end
    return string("Unknown DynEntry Type (0x", string(d_type, base=16), ")")
end

struct ELFDynEntryRef{H <: ELFHandle}
    section_ref::SectionRef{H}
    entry::ELFDynEntry{H}
end

Section(d::ELFDynEntryRef) = d.section_ref
handle(d::ELFDynEntryRef) = handle(Section(d))
deref(d::ELFDynEntryRef) = d.entry
@derefmethod dyn_entry_type(d::ELFDynEntryRef)
@derefmethod dyn_entry_is_string(d::ELFDynEntryRef)
strtab_lookup(d::ELFDynEntryRef) = strtab_lookup(StrTab(d), deref(d).d_val)

function show(io::IO, d::DT) where {DT <: ELFDynEntry}
    print(io, "ELFDynEntry", DT <: ELFDynEntry64 ? " (64 bit)" : "")
    print(io, " $(dyn_entry_type_string(d))")
end

function show(io::IO, d::ELFDynEntryRef)
    show(io, deref(d))
    if dyn_entry_is_string(d)
        print(io, ", \"$(strtab_lookup(d))\"")
    end
end

"""
    StrTab(d::ELFDynEntryRef)

Given an `ELFDynEntryRef`, construct an `ELFStrTab` that refers to the
associated dynamic string table.
"""
function StrTab(d::ELFDynEntryRef{H}) where {H <: ObjectHandle}
    # The `sh_link` field in a ".dynamic" section always points to the
    # dynamic string table (usually called ".dynstr")
    section_idx = deref(Section(d)).sh_link+1
    strtab_section = Sections(handle(d))[section_idx]
    return ELFStrTab{H}(strtab_section)
end

"""
    ELFDynEntries(oh::ELFHandle)

Read all `ELFDynEntry` objects from an ELF object, returning them as an array.
"""
function ELFDynEntries(oh::ELFHandle)
    dyn_section = findfirst(Sections(oh), ".dynamic")

    # Figure out what kind of DynEntry objects we're going to read in
    DT = dyn_entry_type(oh)

    # Seek to the proper spot, and read 'em in until we reach a DT_NULL
    seekstart(dyn_section)
    dts = [ELFDynEntryRef(dyn_section, unpack(oh, DT))]

    while dyn_entry_type(dts[end]) != DT_NULL
        push!(dts, ELFDynEntryRef(dyn_section, unpack(oh, DT)))
    end

    return dts
end

"""
    ELFDynEntries(oh::ELFHandle, kinds::Vector)

Read all `ELFDynEntry` objects from an ELF object, returning them as an array
if they are one of the `kinds` passed in, such as `DT_NEEDED`.
"""
function ELFDynEntries(oh::ELFHandle, kinds::Vector)
    return [d for d in ELFDynEntries(oh) if dyn_entry_type(d) in kinds]
end


text:
# DESCRIPTION
# Model order reduction based on residualization. Keep only rigid-body states.
#
# ------------------------------------------------------------------------------

function ModRed8(A,B)
    idxss=collect(1:8)
    idxfs=[collect(13:31); collect(33:23)]
    Ared=A[idxss,idxss]-A[idxss,idxfs]/A[idxfs,idxfs]*A[idxfs,idxss]
    Bred=B[idxss,1:4]-A[idxss,idxfs]/A[idxfs,idxfs]*B[idxfs,1:4]
    return Ared, Bred
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
__precompile__()

module BeeswaxEsolang

include("beeswax.jl")
include("bwinstructions.jl")
include("bwdebug.jl")
include("prettyprint.jl")

export beeswax
export prettyprint

end     #module


text:
import PDMats.dim
import ConjugatePriors.NormalWishart

immutable MultivariateStudent
    nu::Float64
    mu::Vector{Float64}
    sigma::AbstractPDMat
end

function MultivariateStudent(nu::Float64, mu::Vector{Float64}, sigma::Matrix{Float64})
    return MultivariateStudent(nu, mu, PDMat(sigma))
end

dim(t::MultivariateStudent) = dim(t.sigma)

function lognorm(t::MultivariateStudent)
    return lgamma(t.nu / 2) + (dim(t)/2) * (log(t.nu) + log(pi)) + 0.5 * logdet(t.sigma) - lgamma((t.nu + dim(t)) / 2)
end

function logpdf{Tf <: FloatingPoint}(t::MultivariateStudent, x::Vector{Tf})
    return - (t.nu + dim(t)) / 2 * log(1 + invquad(t.sigma, x - t.mu) / t.nu) - lognorm(t)
end

function pdf{Tf <: FloatingPoint}(t::MultivariateStudent, x::Vector{Tf})
    return exp(logpdf(t, x))
end

function predictive(nw::NormalWishart)
    nu = nw.nu - nw.dim + 1
    return MultivariateStudent(nu, nw.mu, inv(nw.Tchol) * (nw.kappa + 1)/(nw.kappa * nu))
end


text:
using CUDA
using Folds
using FoldsCUDA
using StaticArrays
using Transducers
using Transducers: whencompletebasecase

# TODO: pretty FLoops syntax
function countints_svector_functional(indices, ::Val{n}; ex = PreferParallel()) where {n}

    function init()
        # On initialization, create a `MVector` as a basecase-local
        # sub-histogram buffer `b`:
        zero(MVector{n,Int})
    end

    function inc!(b, i)
        @inbounds b[max(begin, min(i, end))] += 1
        b
    end

    function completebasecase(b)
        # After basecase computing is done, convert the buffer `b` to an
        # immutable value `SVector` to share the value across threads:
        SVector(b)
    end

    function combine(h, b)
        # Cross-thread reduction is simply point-wise addition:
        h .+ b
    end

    rf =
        inc! |>
        wheninit(init) |>
        whencompletebasecase(completebasecase) |>
        whencombine(combine)

    Folds.reduce(rf, indices, ex)
end