text:
using Missings

struct LinkBranch{RT, NL} <: AbstractBranch{RT, NL}
    inout::Tuple{AbstractNode{RT, NL}, AbstractNode{RT, NL}}
    length::Float64
    data::Dict{String, Any}
end

mutable struct LinkNode{RT, NL} <: AbstractNode{RT, NL}
    name::NL
    inbound::Union{AbstractBranch{RT, NL}, Missing}
    other::Vector{AbstractBranch{RT, NL}}
    data::Dict{String, Any}
end

mutable struct LinkTree{RT, NL, N <: LinkNode{RT, NL},
                        B <: LinkBranch{RT, NL}, TD} <:
               AbstractTree{OneTree, RT, NL, N, B}
    name::Union{String, Missing}
    nodes::Dict{NL, N}
    roots::Vector{N}
    branches::Vector{Union{B, Missing}}
    data::Dict{String, Any}
    tipdata::TD
    rootheight::Float64
    isvalid::Union{Bool, Missing}

    function LinkTree{RT, NL, N, B, TD}(tipnames::Vector{NL} = NL[];
        treename::Union{String, Missing} = missing,
        tipdata::TD = TD(), rootheight = NaN) where {RT, NL, N, B, TD}
        tree = new{RT, NL, N, B, TD}(treename, Dict{NL, N}(), Vector{N}(),
                                     Vector{B}(), Dict{String, Any}(),
                                     tipdata, rootheight, missing)
        if !isempty(tipnames)
            createnodes!(tree, tipnames)
        elseif !isempty(tipdata)
            tips = unique(keys(tipdata))
            createnodes!(tree, tips)
        end
        return tree
    end
end
const LB{RT} = LinkBranch{RT, String}
const LN{RT} = LinkNode{RT, String}
const LT{RT, TD} = LinkTree{RT, String, LN{RT}, LB{RT}, TD}
const RootedTree = LT{OneRoot, Dict{String, Any}}
const ManyRootTree = LT{ManyRoots, Dict{String, Any}}
const UnrootedTree = LT{Unrooted, Dict{String, Any}}

# LinkBranch methods
function LinkBranch(from::LinkNode{RT, NL},
                    to::LinkNode{RT, NL},
                    length::Float64 = NaN,
                    data::Dict{String, Any} = Dict{String, Any}()) where
    {RT, NL}
    return LinkBranch{RT, NL}((from, to), length, data)
end

import Phylo.API._src
Phylo.API._src(branch::LinkBranch{RT, NL}) where {RT, NL} =
    branch.inout[1]

import Phylo.API._dst
Phylo.API._dst(branch::LinkBranch{RT, NL}) where {RT, NL} =
    branch.inout[2]

import Phylo.API._getlength
Phylo.API._getlength(branch::LinkBranch{RT, NL}) where {RT, NL} =
    branch.length

import Phylo.API._branchrecordtype
Phylo.API._branchrecordtype(::Type{LinkBranch{RT, NL}}) where {RT, NL} =
    Dict{String, Any}

import Phylo.API._getbranchinfo
Phylo.API._getbranchinfo(branch::LinkBranch{RT, NL}) where {RT, NL} =
    branch.data

import Phylo.API._setbranchinfo!
Phylo.API._setbranchinfo!(branch::LinkBranch{RT, NL},
                          data::Dict{String, Any}) where {RT, NL} =
    branch.data = data

# LinkNode methods
function LinkNode(tree::T, name::NL,
                  data::Dict{String, Any} = Dict{String, Any}()) where
    {RT, NL, N <: LinkNode, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    return N(name, missing, Vector{B}(), data)
end

import Phylo.API._noderecordtype
Phylo.API._noderecordtype(::Type{LinkNode{RT, NL}}) where {RT, NL} =
    Dict{String, Any}

import Phylo.API._getnoderecord
Phylo.API._getnoderecord(node::LinkNode{RT, NL}) where {RT, NL} =
    node.data

import Phylo.API._setnoderecord!
Phylo.API._setnoderecord!(node::LinkNode{RT, NL},
                        data::Dict{String, Any}) where {RT, NL} =
    node.data = data

import Phylo.API._hasinbound
_hasinbound(node::LinkNode{RT, NL}) where {RT <: Rooted, NL} =
    !ismissing(node.inbound)

import Phylo.API._degree
_degree(node::LinkNode{Unrooted, NL}) where {NL} =
    length(node.other)

import Phylo.API._getinbound
_getinbound(node::LinkNode) = node.inbound

import Phylo.API._addinbound!
function _addinbound!(node::LinkNode{RT, NL}, branch::B) where
         {RT <: Rooted, NL, B <: AbstractBranch{RT, NL}}
    !_hasinbound(node) ||
        error("LinkNode already has an inbound connection")
    node.inbound = branch
end

import Phylo.API._removeinbound!
function _removeinbound!(node::LinkNode{RT, NL}, branch::B) where
    {RT <: Rooted, NL, B <: LinkBranch{RT, NL}}
    _hasinbound(node) ||
        error("Node has no inbound connection")
    node.inbound === branch ||
        error("Node has no inbound connection from branch $inbound")
    node.inbound = missing
end

import Phylo.API._getoutbounds
_getoutbounds(node::LinkNode{RT, NL}) where
    {RT <: Rooted, NL} = node.other

import Phylo.API._addoutbound!
function _addoutbound!(node::LinkNode{RT, NL}, branch::B) where
         {RT <: Rooted, NL, B <: LinkBranch{RT, NL}}
    push!(node.other, branch)
end

import Phylo.API._removeoutbound!
function _removeoutbound!(node::LinkNode{RT, NL}, branch::B) where
         {RT <: Rooted, NL, B <: LinkBranch{RT, NL}}
    branch ∈ _getoutbounds(node) ? filter!(p -> p !== branch, node.other) :
         error("Node does not have outbound connection to branch $branch")
end

import Phylo.API._getconnections
_getconnections(node::LinkNode{Unrooted, NL}) where
    {NL} = node.other

import Phylo.API._addconnection!
function _addconnection!(node::LinkNode{Unrooted, NL}, branch::B) where
    {NL, B <: LinkBranch{Unrooted, NL}}
    push!(node.other, branch)
end

import Phylo.API._removeconnection!
function _removeconnection!(node::LinkNode{Unrooted, NL}, branch::B) where
    {NL, B <: LinkBranch{Unrooted, NL}}
    outbound ∈ _getoutbounds(node) ? filter!(p -> p !== branch, node.other) :
         error("Node does not have a connection to branch $branch")
end

# LinkTree methods
const TREENAME = "Tree"
import Phylo.API._treenametype
_treenametype(::Type{T}) where
    {RT, NL, N, B, T <: LinkTree{RT, NL, N, B}} = typeof(TREENAME)
import Phylo.API._gettreename
_gettreename(tree::LinkTree{RT, NL, N, B, TD}) where
    {RT, NL, N, B, TD} = ismissing(tree.name) ? TREENAME : tree.name

import Phylo.API._addnode!
function Phylo.API._addnode!(tree::LinkTree{RT, NL, N, B, TD},
    node::N, name::NL) where {RT, NL, N, B, TD}
    tree.nodes[name] = node
    push!(tree.roots, node)
    isvalid = missing
    return name
end

import Phylo.API._removenode!
function Phylo.API._removenode!(tree::LinkTree{RT, NL, N, B, TD},
                                node::N) where {RT, NL, N, B, TD}
    delete!(tree.nodes, node.name)
    filter!(n -> n !== node, tree.roots)
    isvalid = missing
end

import Phylo.API._deletenode!
_deletenode!(tree::LinkTree{RT, NL, N, B, TD}, name::NL) where
    {RT, NL, N, B, TD} = _deletenode!(tree, tree.nodes[name])
function _deletenode!(tree::LinkTree{RT, NL, N, B, TD}, node::N) where
    {RT <: Rooted, NL, N, B, TD}
    _hasinbound(node) && _deletebranch!(tree, node.inbound)
    while _outdegree(node) > 0
        _deletebranch!(tree, first(node.other))
    end
    _removenode!(tree, node)
end
function _deletenode!(tree::LinkTree{Unrooted, NL, N, B, TD}, branch::B) where
    {NL, N, B, TD}
    while _degree(node) > 0
        _deletebranch!(tree, first(node.other))
    end
    _removenode!(tree, node)
end

import Phylo.API._getroots
_getroots(tree::LinkTree{RT, NL, N, B, TD}) where
    {RT <: Rooted, NL, N, B, TD} = tree.roots
_getroots(tree::LinkTree{Unrooted, NL, N, B, TD}) where
    {NL, N, B, TD} = error("Unrooted trees do not have roots")

import Phylo.API._getnodes
_getnodes(tree::LinkTree{RT, NL, N, B, TD}) where
    {RT, NL, N, B, TD} = values(tree.nodes)

import Phylo.API._getnodenames
_getnodenames(tree::LinkTree{RT, NL, N, B, TD}) where
    {RT, NL, N, B, TD} = keys(tree.nodes)

import Phylo.API._getnodename
_getnodename(::LinkTree{RT, NL, N, B, TD}, node::N) where
    {RT, NL, N, B, TD} = node.name

import Phylo.API._hasnode
_hasnode(tree::LinkTree{RT, NL, N, B, TD}, name::NL) where
    {RT, NL, N, B, TD} = haskey(tree.nodes, name)
_hasnode(tree::LinkTree{RT, NL, N, B, TD}, node::N) where
    {RT, NL, N, B, TD} = node ∈ values(tree.nodes)

import Phylo.API._getnode
_getnode(tree::LinkTree{RT, NL, N, B, TD}, name::NL) where
    {RT, NL, N, B, TD} = tree.nodes[name]

import Phylo.API._getbranches
_getbranches(tree::LinkTree{RT, NL, N, B, TD}) where
    {RT, NL, N, B, TD} = skipmissing(tree.branches)

import Phylo.API._getbranchnames
_getbranchnames(tree::LinkTree{RT, NL, N, B, TD}) where
    {RT, NL, N, B, TD} =
    filter(i -> !ismissing(tree.branches[i]), 1:length(tree.branches))

import Phylo.API._getbranchname
Phylo.API._getbranchname(tree::LinkTree{RT, NL, N, B, TD}, id::Int) where
    {RT, NL, N, B, TD} = id
Phylo.API._getbranchname(tree::LinkTree{RT, NL, N, B, TD}, branch::B) where
    {RT, NL, N, B, TD} = first(indexin([branch], tree.branches))

import Phylo.API._hasbranch
_hasbranch(tree::LinkTree{RT, NL, N, B, TD}, id) where
    {RT, NL, N, B, TD} =
    1 ≤ id ≤ length(tree.branches) && !ismissing(tree.branches[id])
_hasbranch(tree::LinkTree{RT, NL, N, B, TD}, branch::B) where
    {RT, NL, N, B, TD} = branch ∈ tree.branches

import Phylo.API._getbranch
_getbranch(tree::LinkTree{RT, NL, N, B, TD}, id::Int) where
    {RT, NL, N, B, TD} = tree.branches[id]

import Phylo.API._createbranch!
function _createbranch!(tree::LinkTree{RT, NL, N, B, TD},
                        from::LinkNode{RT, NL},
                        to::LinkNode{RT, NL},
                        length::Float64 = NaN,
                        data::Dict{String, Any} = Dict{String, Any}()) where
    {RT <: Rooted, NL, N, B, TD}
    branch = LinkBranch(from, to, length, data)
    _addoutbound!(from, branch)
    _addinbound!(to, branch)
    return _addbranch!(tree, branch)
end
function _createbranch!(tree::LinkTree{Unrooted, NL, N, B, TD},
                        from::LinkNode{Unrooted, NL},
                        to::LinkNode{Unrooted, NL},
                        length::Float64 = NaN,
                        data::Dict{String, Any} = Dict{String, Any}()) where
    {NL, N, B, TD}
    branch = LinkBranch(tree, from, to, length, data)
    _addconnection!(from, branch)
    _addconnection!(to, branch)
    return _addbranch!(tree, branch)
end

import Phylo.API._addbranch!
function _addbranch!(tree::LinkTree{RT, NL, N, B, TD},
                     branch::B) where {RT <: Rooted, NL, N, B, TD}
    id = _newbranchlabel(tree)
    len = length(tree.branches)
    if id <= len
        tree.branches[id] = branch
    elseif id == len + 1
        push!(tree.branches, branch)
    else
        error("Adding branch beyond end of branch vector")
    end

    filter!(n -> n !== _dst(branch), tree.roots)
    isvalid = missing
    return id
end
function _addbranch!(tree::LinkTree{Unrooted, NL, N, B, TD},
                     branch::B, id::Int) where {NL, N, B, TD}
    id = _newbranchlabel(tree)
    tree.branches[id] = branch
    isvalid = missing
    return id
end

import Phylo.API._removebranch!
function Phylo.API._removebranch!(tree::LinkTree{RT, NL, N, B, TD},
                                  branch::B,
                                  id::Int = first(indexin([branch],
                                                  tree.branches))) where
    {RT <: Rooted, NL, N, B, TD}
    tree.branches[id] = missing
    push!(tree.roots, _dst(branch))
    isvalid = missing
end
function Phylo.API._removebranch!(tree::LinkTree{Unrooted, NL, N, B, TD},
                                  branch::B,
                                  id::Int = first(indexin([branch],
                                                  tree.branches))) where
    {NL, N, B, TD}
    tree.branches[id] = missing
    isvalid = missing
end

import Phylo.API._deletebranch!
_deletebranch!(tree::LinkTree{RT, NL, N, B, TD}, id::Int) where
    {RT, NL, N, B, TD} = _deletebranch!(tree, tree.branches[id], id)
function _deletebranch!(tree::LinkTree{RT, NL, N, B, TD}, branch::B,
                        id::Int = first(indexin([branch], tree.branches))) where
    {RT <: Rooted, NL, N, B, TD}
    _removeoutbound!(branch.inout[1], branch)
    _removeinbound!(branch.inout[2], branch)
    _removebranch!(tree, branch, id)
end
function _deletebranch!(tree::LinkTree{Unrooted, NL, N, B, TD}, branch::B,
                        id::Int = first(indexin([branch], tree.branches))) where
    {NL, N, B, TD}
    _removeconnection!(branch.inout[1], branch)
    _removeconnection!(branch.inout[2], branch)
    _removebranch!(tree, branch, id)
end

import Phylo.API._createnode!
function _createnode!(tree::LinkTree{RT, NL, N, B, TD}, nodename::NL,
                      data::Dict{String, Any} = Dict{String, Any}()) where
    {RT <: Rooted, NL, N, B, TD}
    node = LinkNode(tree, nodename, data)
    return _addnode!(tree, node, nodename)
end

import Phylo.API._getleafinfo
_getleafinfo(tree::LinkTree{RT, NL, N, B, TD}) where
    {RT, NL, N, B, TD} = tree.tipdata

import Phylo.API._noderecordtype
Phylo.API._noderecordtype(::Type{LinkTree{RT, NL, N, B, TD}}) where
    {RT, NL, N, B, TD} = _noderecordtype(N)

import Phylo.API._getnoderecord
Phylo.API._getnoderecord(tree::LinkTree{RT, NL, N, B, TD}, name::NL) where
    {RT, NL, N, B, TD} = _getnoderecord(_getnode(tree, name))
Phylo.API._getnoderecord(tree::LinkTree{RT, NL, N, B, TD}, node::N) where
    {RT, NL, N, B, TD} = _getnoderecord(node)

import Phylo.API._branchrecordtype
Phylo.API._branchrecordtype(::Type{LinkTree{RT, NL, N, B, TD}}) where
    {RT, NL, N, B, TD} = _branchrecordtype(N)

import Phylo.API._getbranchinfo
Phylo.API._getbranchinfo(tree::LinkTree{RT, NL, N, B, TD}, id::Int) where
    {RT, NL, N, B, TD} = _getbranchinfo(_getbranch(tree, id))
Phylo.API._getbranchinfo(tree::LinkTree{RT, NL, N, B, TD}, branch::B) where
    {RT, NL, N, B, TD} = _getbranchinfo(branch)


text:
module MapWindow

using DataStructures, TiledIteration
using ..ImageFiltering: BorderSpecAny, Pad, Fill, Inner,
    borderinstance, _interior, padindex, imfilter
using Base: Indices, tail
using Statistics

export mapwindow, mapwindow!

"""
    mapwindow(f, img, window; [border="replicate"], [indices=axes(img)]) -> imgf

Apply `f` to sliding windows of `img`, with window size or axes
specified by `window`. For example, `mapwindow(median!, img, window)`
returns an `Array` of values similar to `img` (median-filtered, of
course), whereas `mapwindow(extrema, img, window)` returns an `Array`
of `(min,max)` tuples over a window of size `window` centered on each
point of `img`.

The function `f` receives a buffer `buf` for the window of data
surrounding the current point. If `window` is specified as a
Dims-tuple (tuple-of-integers), then all the integers must be odd and
the window is centered around the current image point. For example, if
`window=(3,3)`, then `f` will receive an Array `buf` corresponding to
offsets `(-1:1, -1:1)` from the `imgf[i,j]` for which this is
currently being computed. Alternatively, `window` can be a tuple of
AbstractUnitRanges, in which case the specified ranges are used for
`buf`; this allows you to use asymmetric windows if needed.

`border` specifies how the edges of `img` should be handled; see
`imfilter` for details.

Finally `indices` allows to omit unnecessary computations, if you want to do things
like `mapwindow` on a subimage, or a strided variant of mapwindow.
It works as follows:
```julia
mapwindow(f, img, window, indices=(2:5, 1:2:7)) == mapwindow(f,img,window)[2:5, 1:2:7]
```
Except more efficiently because it omits computation of the unused values.

Because the data in the buffer `buf` that is received by `f` is copied from `img`, and the
buffer's memory is reused, `f` should not return references to `buf`.
This
```julia
f = buf->copy(buf) # as opposed to f = buf->buf
mapwindow(f, img, window, indices=(2:5, 1:2:7))
```
would work as expected.

For functions that can only take `AbstractVector` inputs, you might have to
first specialize `default_shape`:

```julia
f = v->quantile(v, 0.75)
ImageFiltering.MapWindow.default_shape(::typeof(f)) = vec
```

and then `mapwindow(f, img, (m,n))` should filter at the 75th quantile.

See also: [`imfilter`](@ref).
"""
function mapwindow(f, img, window; border="replicate",
                   indices=default_imginds(img, window, border), callmode=:copy!)
    if callmode != :copy!
        error("Only callmode=:copy! is currently supported")
    end
    _mapwindow(replace_function(f),
              img,
              resolve_window(window),
              resolve_border(border),
              resolve_imginds(indices))
end

function default_imginds(img, window, border)
    axes(img)
end
function default_imginds(img, window, border::Inner)
    imginds = axes(img)
    win = resolve_window(window)
    indind = _indices_of_interiour_indices(imginds, imginds, win)
    map((II, r) -> r[II], indind, imginds)
end

function _mapwindow(f, img, window, border, imginds)
    out = allocate_output(f, img, window, border, imginds)
    mapwindow_kernel!(f, out, img, window, border, imginds)
end

"""
    mapwindow!(f, out, img, window; border="replicate", indices=axes(img))

Variant of [mapwindow](@ref), with preallocated output.
If `out` and `img` have overlapping memory regions, behaviour is undefined.
"""
function mapwindow!(f, out, img, window; border="replicate",
                    indices=default_imginds(img, window, border))
    mapwindow_kernel!(replace_function(f),
              out,
              img,
              resolve_window(window),
              resolve_border(border),
              resolve_imginds(indices))
end

function median_fast!(v)
    # median! calls partialsort! which has keyword arguments. Keyword arguments are slow.
    # This replaces median! with a more efficient implementation free of keyword arguments.
    inds = axes(v,1)
    Statistics.middle(Base.partialsort!(v, (first(inds)+last(inds))÷2, Base.Order.ForwardOrdering()))
end

replace_function(f) = f
replace_function(::typeof(median!)) = median_fast!

function resolve_window(window::Dims)
    all(isodd(w) for w in window) || error("entries in window must be odd, got $window")
    halfsize = map(w->w>>1, window)
    map(h -> -h:h,halfsize)
end
function resolve_window(window::Integer)
    isodd(window) || error("window must be odd, got $window")
    h = window>>1
    (-h:h,)
end
resolve_window(window::AbstractArray) = resolve_window((window...,))
resolve_window(window::AbstractUnitRange) = (window,)
resolve_window(window::Indices) = window

resolve_border(border::AbstractString) = borderinstance(border)
resolve_border(border::BorderSpecAny) = border

resolve_imginds(r::AbstractRange) = (r,)
resolve_imginds(imginds) = imginds

abstract type _IndexTransformer end

struct _AffineTransformer{N} <: _IndexTransformer
    offset::NTuple{N,Int}
    stride::NTuple{N,Int}
end
@inline function Base.getindex(t::_AffineTransformer, inds::CartesianIndex)
    CartesianIndex(t.offset .+ t.stride .* inds.I)
end

struct _OffsetTransformer{N} <: _IndexTransformer
    offset::NTuple{N,Int}
end
@inline function Base.getindex(t::_OffsetTransformer, inds::CartesianIndex)
    CartesianIndex(t.offset .+ inds.I)
end

struct _IdentityTransformer <: _IndexTransformer end
@inline Base.getindex(t::_IdentityTransformer, inds) = inds

function _IndexTransformer(ranges)
    stride = map(step, ranges)
    offset1 = map(first, ranges)
    offset = offset1 .- stride
    _AffineTransformer(offset, stride)
end

function _IndexTransformer(ranges::NTuple{N,Base.OneTo}) where {N}
    _IdentityTransformer()
end

function _IndexTransformer(ranges::NTuple{N,AbstractUnitRange}) where {N}
    offset1 = map(first, ranges)
    offset = offset1 .- 1
    _OffsetTransformer(offset)
end

compute_output_range(r::AbstractUnitRange) = r
compute_output_range(r::AbstractRange) = Base.OneTo(length(r))

function compute_output_indices(imginds)
    ranges = map(compute_output_range, imginds)
    # Base.similar does not like if some but not all ranges are Base.OneTo
    homogenize(ranges)
end
homogenize(ranges::NTuple{N, AbstractRange}) where {N}   = map(r-> first(r):step(r):last(r), ranges)
homogenize(ranges::NTuple{N, AbstractUnitRange}) where{N} = map(r-> first(r):last(r), ranges)
homogenize(ranges::NTuple{N, Base.OneTo}) where {N} = ranges

# Return indices of elements of `r` that are also elements of `full`.
function _intersectionindices(full::AbstractUnitRange, r::AbstractRange)
    r_sub = intersect(full, r)
    if isempty(r_sub)
        ret = 1:0
    else
        ret = _indexof(r,first(r_sub)):_indexof(r,last(r_sub))
    end
    @assert intersect(full, r) == r[ret]
    ret
end

function _indexof(r::AbstractRange, x)
    T = eltype(r)
    @assert x ∈ r
    i = one(T) + T((x - first(r)) / step(r))
    @assert r[i] == x
    i
end

function _indices_of_interiour_range(
        fullimgr::AbstractUnitRange,
        imgr::AbstractRange,
        kerr::AbstractRange)
    kmin, kmax = extrema(kerr)
    idx1 = _intersectionindices(fullimgr, kmin .+ imgr)
    idx2 = _intersectionindices(fullimgr, kmax .+ imgr)
    idx = intersect(idx1, idx2)
    @assert imgr[idx] .+ kmin ⊆ fullimgr
    @assert imgr[idx] .+ kmax ⊆ fullimgr
    idx
end

function _indices_of_interiour_indices(fullimginds, imginds, kerinds)
    map(_indices_of_interiour_range, fullimginds, imginds, kerinds)
end

function allocate_output(f, img, window, border, imginds)
    T = compute_output_eltype(f, img, window, border, imginds)
    outinds = compute_output_indices(imginds)
    similar(img, T, outinds)
end

function allocate_buffer(f, img, window)
    T = eltype(img)
    buf = Array{T}(undef,map(length, window))
    bufrs = default_shape(f)(buf)
    buf, bufrs
end

function compute_output_eltype(f, img, window, border, imginds)
    buf, bufrs = allocate_buffer(f, img, window)
    make_buffer_values_realistic!(buf, img, window, border, imginds)
    typeof(f(bufrs))
end

function make_buffer_values_realistic!(buf, img, window, border::Inner, imginds)
    x = oneunit(eltype(img))
    fill!(buf, x)
end

function make_buffer_values_realistic!(buf, img, window, border, imginds)
    Iimg = CartesianIndex(map(first, imginds))
    offset = CartesianIndex(map(w->first(w)-1, window))
    copy_win!(buf, img, Iimg, border, offset)
end

function mapwindow_kernel!(f,
                    out::AbstractArray{S,N},
                    img::AbstractArray{T,N},
                    window::NTuple{N,AbstractUnitRange},
                    border::BorderSpecAny,
                    imginds::NTuple{N, AbstractRange}) where {S,T,N}

    @assert map(length, imginds) == map(length, axes(out))

    indind_full = map(r -> Base.OneTo(length(r)), imginds)
    indind_inner = _indices_of_interiour_indices(axes(img), imginds, window)
    Rindind_full = CartesianIndices(indind_full)
    Rindind_inner = CartesianIndices(indind_inner)

    outindtrafo = _IndexTransformer(axes(out))
    imgindtrafo = _IndexTransformer(imginds)

    buf, bufrs = allocate_buffer(f, img, window)
    Rbuf = CartesianIndices(size(buf))
    for II ∈ Rindind_inner
        Iimg = imgindtrafo[II]
        Iout = outindtrafo[II]
        Rwin = CartesianIndices(map((w,o) -> w .+ o, window, Tuple(Iimg)))
        copyto!(buf, Rbuf, img, Rwin)
        out[Iout] = f(bufrs)
    end
    # Now pick up the edge points we skipped over above
    Rindind_edge = EdgeIterator(Rindind_full, Rindind_inner)
    offset = CartesianIndex(map(w->first(w)-1, window))
    for II ∈ Rindind_edge
        Iimg = imgindtrafo[II]
        Iout = outindtrafo[II]
        copy_win!(buf, img, Iimg, border, offset)
        out[Iout] = f(bufrs)
    end
    out
end

# For copying along the edge of the image
function copy_win!(buf::AbstractArray, img, I, border::Pad, offset)
    win_inds = map((x,y)->x .+ y, axes(buf), Tuple(I) .+ Tuple(offset))
    win_img_inds = map(intersect, axes(img), win_inds)
    padinds = map((inner,outer)->padindex(border, inner, outer), win_img_inds, win_inds)
    docopy!(buf, img, padinds)
    buf
end
docopy!(buf, img, padinds::NTuple{1}) = buf[:] = view(img, padinds[1])
docopy!(buf, img, padinds::NTuple{2}) = buf[:,:] = view(img, padinds[1], padinds[2])
docopy!(buf, img, padinds::NTuple{3}) = buf[:,:,:] = view(img, padinds[1], padinds[2], padinds[3])
@inline function docopy!(buf, img, padinds::NTuple{N}) where N
    colons = ntuple(d->Colon(), Val{N}())
    buf[colons...] = view(img, padinds...)
end

function copy_win!(buf::AbstractArray, img, I, border::Fill, offset)
    R = CartesianIndices(axes(img))
    Ioff = I+offset
    for J in CartesianIndices(axes(buf))
        K = Ioff+J
        buf[J] = K ∈ R ? img[K] : convert(eltype(img), border.value)
    end
    buf
end

### Optimizations for particular window-functions

mapwindow(::typeof(extrema), A::AbstractArray, window::Dims) = extrema_filter(A, window)
mapwindow(::typeof(extrema), A::AbstractVector, window::Integer) = extrema_filter(A, (window,))

# Max-min filter

# This is an implementation of the Lemire max-min filter
# http://arxiv.org/abs/cs.DS/0610046

# Monotonic wedge
struct Wedge{T}
    L::CircularDeque{T}
    U::CircularDeque{T}
end
Wedge{T}(n::Integer) where {T} = Wedge(CircularDeque{T}(n), CircularDeque{T}(n))

function Base.push!(W::Wedge, i::Integer)
    push!(W.L, i)
    push!(W.U, i)
    W
end

function addtoback!(W::Wedge, A, i, J)
    mn, mx = A[i, J]
    @inbounds while !isempty(W.L) && mn < A[back(W.L), J][1]
        pop!(W.L)
    end
    @inbounds while !isempty(W.U) && mx > A[back(W.U), J][2]
        pop!(W.U)
    end
    push!(W.L, i)
    push!(W.U, i)
    W
end

function Base.empty!(W::Wedge)
    empty!(W.L)
    empty!(W.U)
    W
end

@inline function getextrema(A, W::Wedge, J)
    (A[front(W.L), J][1], A[front(W.U), J][2])
end

"""
    extrema_filter(A, window) --> Array{(min,max)}

Calculate the running min/max over a window of width `window[d]` along
dimension `d`, centered on the current point. The returned array has
the same axes as the input `A`.
"""
function extrema_filter(A::AbstractArray{T,N}, window::NTuple{N,Integer}) where {T,N}
    _extrema_filter!([(a,a) for a in A], window...)
end
extrema_filter(A::AbstractArray, window::AbstractArray) = extrema_filter(A, (window...,))
extrema_filter(A::AbstractArray, window) = error("`window` must have the same number of entries as dimensions of `A`")

extrema_filter(A::AbstractArray{T,N}, window::Integer) where {T,N} = extrema_filter(A, ntuple(d->window, Val{N}))

function _extrema_filter!(A::Array, w1, w...)
    if w1 > 1
        a = first(A)
        if w1 <= 20
            cache = ntuple(i->a, w1>>1)
            _extrema_filter1!(A, w1, cache)
        else
            n = w1>>1
            cache = CircularDeque{typeof(a)}(n)
            for i = 1:n
                push!(cache, a)
            end
            _extrema_filter1!(A, w1, cache)
        end
    end
    if ndims(A) > 1
        _extrema_filter!(permutedims(A, [2:ndims(A);1]), w...)
    else
        return A
    end
end
_extrema_filter!(A::Array) = A

# Extrema-filtering along "columns" (dimension 1). This implements Lemire
# Algorithm 1, with the following modifications:
#   - multidimensional array support by looping over trailing dimensions
#   - working with min/max pairs rather than plain values, to
#     facilitate multidimensional processing
#   - output for all points of the array, handling the edges as max-min
#     over halfwindow on either side
function _extrema_filter1!(A::AbstractArray{Tuple{T,T}}, window::Int, cache) where T
    # Initialise the internal wedges
    # U[1], L[1] are the location of the global (within the window) maximum and minimum
    # U[2], L[2] are the maximum and minimum over (U1, end] and (L1, end], respectively
    W = Wedge{Int}(window+1)
    tmp = Array{Tuple{T,T}}(undef, window)
    c = z = first(cache)

    inds = axes(A)
    inds1 = inds[1]
    halfwindow = window>>1
    iw = min(last(inds1), first(inds1)+window-1)
    for J in CartesianIndices(tail(inds))
        empty!(W)
        # Leading edge. We can't overwrite any values yet in A because
        # we'll need them again in later computations.
        for i = first(inds1):iw
            addtoback!(W, A, i, J)
            c, cache = cyclecache(cache, getextrema(A, W, J))
        end
        # Process the rest of the "column"
        for i = iw+1:last(inds1)
            A[i-window, J] = c
            if i == window+front(W.U)
                popfirst!(W.U)
            end
            if i == window+front(W.L)
                popfirst!(W.L)
            end
            addtoback!(W, A, i, J)
            c, cache = cyclecache(cache, getextrema(A, W, J))
        end
        for i = last(inds1)-window+1:last(inds1)-1
            if i >= first(inds1)
                A[i, J] = c
            end
            if i == front(W.U)
                popfirst!(W.U)
            end
            if i == front(W.L)
                popfirst!(W.L)
            end
            c, cache = cyclecache(cache, getextrema(A, W, J))
        end
        A[last(inds1), J] = c
    end
    A
end

# This is slightly faster than a circular buffer
@inline cyclecache(b::Tuple, x) = b[1], (Base.tail(b)..., x)
@inline function cyclecache(b::CircularDeque, x)
    ret1 = popfirst!(b)
    push!(b, x)
    return ret1, b
end

default_shape(::Any) = identity
default_shape(::typeof(median_fast!)) = vec

## Deprecations
function mapwindow(f, img, window, border)
    Base.depwarn("mapwindow(f, img, window, $border) is deprecated, use `mapwindow(f, img, window, border=$border)` instead.", :mapwindow)
    mapwindow(f,img,window,border=border)
end

function mapwindow(f, img, window, border, indices)
    Base.depwarn("mapwindow(f, img, window, $border, $indices) is deprecated, use `mapwindow(f, img, window, border=$border, indices=$indices)` instead.", :mapwindow)
    mapwindow(f,img,window,border=border,indices=indices)
end

function mapwindow!(f, out, img, window, border)
    Base.depwarn("mapwindow!(f, out, img, window, $border) is deprecated, use `mapwindow!(f, out, img, window, border=$border)` instead.", :mapwindow!)
    mapwindow!(f,out,img,window,border=border)
end

function mapwindow!(f, out, img, window, border, indices)
    Base.depwarn("mapwindow!(f, out, img, window, $border, $indices) is deprecated, use `mapwindow!(f, out, img, window, border=$border, indices=$indices)` instead.", :mapwindow!)
    mapwindow!(f,out,img,window,border=border,indices=indices)
end

end

ext: jl


text:
import ProximalOperators: gradient, gradient!, prox!, prox

struct Zero end

(f::Zero)(x) = 0.0

function gradient!(y, f::Zero, x)
    fill!(y,0.0)
    return 0.0
end

function gradient(f::Zero, x)
    y = zero(x)
    return y, 0.0
end

function prox!(y, f::Zero, x, gamma)
    y .= x
    return 0.0
end

function prox(f::Zero, x, gamma)
    y = copy(x)
    return y, 0.0
end


text:
module M2
using ConcreteInterfaces

abstract type A end

struct B <: A
    x::Int
end

struct C <: A
    x::Float64
end

@interface Foo{T}(self::A) begin
    foo(y::T)::Base.promote_op(+, T, Float64) = self.x + y
    bar(y::T)::T = self.x
end
end

@testset "simple interface with parent" begin
    b = M2.B(1)
    c = M2.C(2.0)
    @test isa(@inferred(M2.Foo{Float64}(b)), M2.Foo{Float64, Float64})
    F = M2.Foo{Float64}.([b, c])
    @test F[1].self === b
    @test F[2].self === c
    @test eltype(F) == M2.Foo{Float64, Float64}
    @test @inferred(M2.foo(F[1], 2.0)) == 3.0
    @test M2.foo.(F, 2.0) == [3.0, 4.0]
    y = zeros(2)
    @test @wrappedallocs(y .= M2.foo.(F, 2.0)) == 0
end


text:
using SparseArrays
const SA = SparseArrays

"""
    τ = freezing_time(x, θ)

computes the hitting time of a 1d particle with
constant velocity `θ` to hit 0 given the position `x`
"""
function freezing_time(x, θ, F::Union{BouncyParticle, ZigZag})
    if θ*x >= 0 # sic!
        return Inf
    else
        return -x/θ
    end
end

"""
    t, x, θ = ssmove_forward!(t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})

moves forward only the non_frozen particles
"""
function ssmove_forward!(t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})
    for i in eachindex(x)
        if θ[i] != 0.0
            t[i], x[i] = t′, x[i] + θ[i]*(t′ - t[i])
        end
    end
    t, x, θ
end
"""
    t, x, θ = ssmove_forward!(G, i, t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})

moves forward only the non_frozen particles neighbours of i
"""
function ssmove_forward!(G, i, t, x, θ, t′, Z::Union{BouncyParticle, ZigZag})
    nhd = neighbours(G, i)
    for i in nhd
        if θ[i] != 0.0
            t[i], x[i] = t′, x[i] + θ[i]*(t′ - t[i])
        end
    end
    t, x, θ
end

"""
    queue_time!(Q, t, x, θ, i, b, f, Z::ZigZag)

Computes the (proposed) reflection time and the freezing time of the
ith coordinate and enqueue the first one. `f[i] = true` if the next
time is a freezing time.
"""
function queue_time!(Q, t, x, θ, i, b, f, Z::ZigZag)
    trefl = poisson_time(b[i], rand())
    tfreeze = freezing_time(x[i], θ[i], Z)
    if tfreeze <= trefl
        f[i] = true
        Q[i] = t[i] + tfreeze
    else
        f[i] = false
        Q[i] = t[i] + trefl
    end
    return Q
end

"""
    sspdmp_inner!(Ξ, G, G1, G2, ∇ϕ, t, x, θ, Q, c, b, t_old, f, θf, (acc, num),
            F::ZigZag, κ, args...; strong_upperbounds = false, factor=1.5, adapt=false)

Inner loop of the sticky ZigZag sampler. `G[i]` are indices which have to be moved, `G1[i]` is the set of indices used to derive the
bounding rate λbar_i and `G2` are the indices k in A_j for all j : i in Aj (neighbours of neighbours)

Is assumed that ∇ϕ[x, i] is function of x_i
with i in G[i] or that ∇ϕ takes care of moving .
"""
function sspdmp_inner!(Ξ, G, G1, G2, ∇ϕ, t, x, θ, Q, c, b, t_old, f, θf, (acc, num),
        F::ZigZag, κ, args...; reversible=false, strong_upperbounds = false, factor=1.5, adapt=false)
    n = length(x)
    # f[i] is true if the next event will be a freeze
    while true
        ii, t′ = peek(Q)
        refresh = ii > n
        i = ii - refresh*n
        refresh && error("refreshment not implemented")
        if f[i] # case 1) to be frozen
            t, x, θ = smove_forward!(i, t, x, θ, t′, F) # move only coordinate i
            if abs(x[i]) > 1e-8
                error("x[i] = $(x[i]) !≈ 0")
            end
            x[i] = -0*θ[i]
            θf[i], θ[i] = θ[i], 0.0 # stop and save speed
            t_old[i] = t[i]
            f[i] = false
            Q[i] = t[i] - log(rand())/κ[i]
            if !strong_upperbounds
                t, x, θ = ssmove_forward!(G, i, t, x, θ, t′, F)
                t, x, θ = ssmove_forward!(G2, i, t, x, θ, t′, F)
                for j in neighbours(G1, i)
                    if θ[j] != 0 # only non-frozen, especially not i
                        b[j] = ab(G1, j, x, θ, c, F, args...)
                        t_old[j] = t[j]
                        Q = queue_time!(Q, t, x, θ, j, b, f, F)
                    end
                end
            end
        elseif x[i] == 0 && θ[i] == 0 # case 2) was frozen
            t[i] = t′ # equivalent to t, x, θ = smove_forward!(i, t, x, θ, t′, F) # move only coordinate i
            θ[i], θf[i] = θf[i], 0.0 # unfreeze, restore speed
            if reversible
                θ[i] *= rand((-1,1))
            end
            t_old[i] = t[i]
            t, x, θ = ssmove_forward!(G, i, t, x, θ, t′, F) # neighbours
            t, x, θ = ssmove_forward!(G2, i, t, x, θ, t′, F) # neighbours of neightbours \ neighbours
            for j in neighbours(G1, i)
                if θ[j] != 0 # only non-frozen, including i # check!
                    b[j] = ab(G1, j, x, θ, c, F, args...)
                    t_old[j] = t[j]
                    Q = queue_time!(Q, t, x, θ, j, b, f, F)
                end
            end
        else # was either a reflection time or an event time from the upper bound
            t, x, θ = ssmove_forward!(G, i, t, x, θ, t′, F) # neighbours
            # do it here so ∇ϕ is right event without self moving
            ∇ϕi = ∇ϕ_(∇ϕ, t, x, θ, i, t′, F, args...)
            l, lb = sλ(∇ϕi, i, x, θ, F), sλ̄(b[i], t[i] - t_old[i])
            num += 1
            if rand()*lb < l # was a reflection time
                acc += 1
                if l > lb
                    !adapt && error("Tuning parameter `c` too small. l/lb = $(l/lb)")
                    acc = num = 0
                    adapt!(c, i, factor)
                end
                # have not moved yet
                t, x, θ = ssmove_forward!(G2, i, t, x, θ, t′, F) # neighbours of neightbours \ neighbours
                θ = reflect!(i, ∇ϕi, x, θ, F)
                for j in neighbours(G1, i)
                    if θ[j] != 0
                        b[j] = ab(G, j, x, θ, c, F, args...)
                        t_old[j] = t[j]
                        queue_time!(Q, t, x, θ, j, b, f, F)
                    end
                end
            else # was an event time from upperbound -> nothing happens
                b[i] = ab(G1, i, x, θ, c, F, args...)
                t_old[i] = t[i]
                queue_time!(Q, t, x, θ, i, b, f, F)
                continue
            end
        end
        push!(Ξ, event(i, t, x, θ, F))
        return t, x, θ, t′, (acc, num), c, b, t_old
    end
end

function sspdmp(∇ϕ, t0, x0, θ0, T, c, G, F::ZigZag, κ, args...; reversible=false,strong_upperbounds = false,
        factor=1.5, adapt=false, progress=false, progress_stops = 20)
    n = length(x0)
    t′ = t0
    t = fill(t′, size(θ0)...)
    t_old = copy(t)
    f = zeros(Bool, n)
    Γ = sparse(F.Γ)
    G1 = [i => rowvals(F.Γ)[nzrange(F.Γ, i)] for i in eachindex(θ0)]
    if G === nothing
        G = G1
    end
    @assert all(a.second ⊇ b.second for (a,b) in zip(G, G1))
    G2 = [i => setdiff(union((G1[j].second for j in G1[i].second)...), G[i].second) for i in eachindex(G1)]
    x, θ = copy(x0), copy(θ0)
    θf = zero(θ)
    num = acc = 0
    Q = SPriorityQueue{Int,Float64}()
    b = [ab(G1, i, x, θ, c, F, args...) for i in eachindex(θ)]
    for i in eachindex(θ)
        trefl = poisson_time(b[i], rand())
        tfreez = freezing_time(x[i], θ[i], F)
        if trefl > tfreez
            f[i] = true
            enqueue!(Q, i => t0 + tfreez)
        else
            f[i] = false
            enqueue!(Q, i => t0 + trefl)
        end
    end
    if hasrefresh(F)
        for i in eachindex(θ)
            enqueue!(Q, (n + i) => waiting_time_ref(F))
        end
    end
    Ξ = Trace(t0, x0, θ0, F)
    if progress
        prg = Progress(progress_stops, 1)
    else
        prg = missing
    end
    stops = ismissing(prg) ? 0 : max(prg.n - 1, 0) # allow one stop for cleanup
    tstop = T/stops
    while t′ < T
        t, x, θ, t′, (acc, num), c,  b, t_old = sspdmp_inner!(Ξ, G, G1, G2, ∇ϕ, t, x, θ, Q,
                    c, b, t_old, f, θf, (acc, num), F, κ, args...; reversible=reversible,
                    strong_upperbounds = strong_upperbounds , factor=factor,
                    adapt=adapt)
        if t′ > tstop
            tstop += T/stops
            next!(prg)
        end
    end
    ismissing(prg) || ProgressMeter.finish!(prg)
    #t, x, θ = ssmove_forward!(t, x, θ, T, F)
    Ξ, (t, x, θ), (acc, num), c
end

sspdmp(∇ϕ, t0, x0, θ0, T, c, F, κ, args...; kwargs...) = sspdmp(∇ϕ, t0, x0, θ0, T, c, nothing, F, κ, args...;  kwargs...)


text:
using ManifoldsBase: DefaultManifold, ℝ

struct PlaneManifold <: AbstractEmbeddedManifold{ℝ,TransparentIsometricEmbedding} end

ManifoldsBase.decorated_manifold(::PlaneManifold) = ManifoldsBase.DefaultManifold(1, 3)
ManifoldsBase.base_manifold(::PlaneManifold) = ManifoldsBase.DefaultManifold(2)

ManifoldsBase.project!(::PlaneManifold, q, p) = (q .= [p[1] p[2] 0.0])
ManifoldsBase.project!(::PlaneManifold, Y, p, X) = (Y .= [X[1] X[2] 0.0])

struct AnotherPlaneManifold <: AbstractEmbeddedManifold{ℝ,DefaultIsometricEmbeddingType} end

ManifoldsBase.decorated_manifold(::AnotherPlaneManifold) = ManifoldsBase.DefaultManifold(3)
ManifoldsBase.base_manifold(::AnotherPlaneManifold) = ManifoldsBase.DefaultManifold(2)

function ManifoldsBase.embed!(::AnotherPlaneManifold, q, p)
    q[1:2] .= p
    q[3] = 0
    return q
end
function ManifoldsBase.embed!(::AnotherPlaneManifold, Y, p, X)
    Y[1:2] .= X
    Y[3] = 0
    return Y
end
function ManifoldsBase.project!(::AnotherPlaneManifold, q, p)
    return q .= [p[1], p[2]]
end
function ManifoldsBase.project!(::AnotherPlaneManifold, Y, p, X)
    return Y .= [X[1], X[2]]
end

function ManifoldsBase.embed!(
    ::EmbeddedManifold{𝔽,DefaultManifold{nL,𝔽},DefaultManifold{mL,𝔽2}},
    q,
    p,
) where {nL,mL,𝔽,𝔽2}
    n = size(p)
    ln = length(n)
    m = size(q)
    lm = length(m)
    (length(n) > length(m)) && throw(
        DomainError(
            "Invalid embedding, since Euclidean dimension ($(n)) is longer than embedding dimension $(m).",
        ),
    )
    any(n .> m[1:ln]) && throw(
        DomainError(
            "Invalid embedding, since Euclidean dimension ($(n)) has entry larger than embedding dimensions ($(m)).",
        ),
    )
    fill!(q, 0)
    q[map(ind_n -> Base.OneTo(ind_n), n)..., ntuple(_ -> 1, lm - ln)...] .= p
    return q
end

function ManifoldsBase.project!(
    ::EmbeddedManifold{𝔽,DefaultManifold{nL,𝔽},DefaultManifold{mL,𝔽2}},
    q,
    p,
) where {nL,mL,𝔽,𝔽2}
    n = size(p)
    ln = length(n)
    m = size(q)
    lm = length(m)
    (length(n) < length(m)) && throw(
        DomainError(
            "Invalid embedding, since Euclidean dimension ($(n)) is longer than embedding dimension $(m).",
        ),
    )
    any(n .< m[1:ln]) && throw(
        DomainError(
            "Invalid embedding, since Euclidean dimension ($(n)) has entry larger than embedding dimensions ($(m)).",
        ),
    )
    #  fill q with the „top left edge“ of p.
    q .= p[map(i -> Base.OneTo(i), m)..., ntuple(_ -> 1, lm - ln)...]
    return q
end

struct NotImplementedEmbeddedManifold <:
       AbstractEmbeddedManifold{ℝ,TransparentIsometricEmbedding} end
function ManifoldsBase.decorated_manifold(::NotImplementedEmbeddedManifold)
    return ManifoldsBase.DefaultManifold(2)
end
function ManifoldsBase.base_manifold(::NotImplementedEmbeddedManifold)
    return ManifoldsBase.DefaultManifold(2)
end

struct NotImplementedEmbeddedManifold2 <:
       AbstractEmbeddedManifold{ℝ,DefaultIsometricEmbeddingType} end

struct NotImplementedEmbeddedManifold3 <: AbstractEmbeddedManifold{ℝ,DefaultEmbeddingType} end

@testset "Embedded Manifolds" begin
    @testset "EmbeddedManifold basic tests" begin
        M = EmbeddedManifold(
            ManifoldsBase.DefaultManifold(2),
            ManifoldsBase.DefaultManifold(3),
        )
        @test repr(M) ==
              "EmbeddedManifold($(sprint(show, M.manifold)), $(sprint(show, M.embedding)))"
        @test base_manifold(M) == ManifoldsBase.DefaultManifold(2)
        @test get_embedding(M) == ManifoldsBase.DefaultManifold(3)
        @test ManifoldsBase.decorated_manifold(M) == ManifoldsBase.DefaultManifold(3)
        @test ManifoldsBase.default_decorator_dispatch(M) === Val(true)
    end

    @testset "PlaneManifold" begin
        M = PlaneManifold()
        @test repr(M) == "PlaneManifold()"
        @test ManifoldsBase.default_decorator_dispatch(M) === Val(false)
        @test ManifoldsBase.default_embedding_dispatch(M) === Val(false)
        @test get_embedding(M) == ManifoldsBase.DefaultManifold(1, 3)
        # Check fallbacks to check embed->check_manifoldpoint Defaults
        @test_throws DomainError is_manifold_point(M, [1, 0, 0], true)
        @test_throws DomainError is_manifold_point(M, [1 0], true)
        @test is_manifold_point(M, [1 0 0], true)
        @test_throws DomainError is_tangent_vector(M, [1 0 0], [1], true)
        @test_throws DomainError is_tangent_vector(M, [1 0 0], [0 0 0 0], true)
        @test is_tangent_vector(M, [1 0 0], [1 0 1], true)
        p = [1.0 1.0 0.0]
        q = [1.0 0.0 0.0]
        X = q - p
        @test check_size(M, p) === nothing
        @test check_size(M, p, X) === nothing
        @test check_size(M, [1, 2]) isa DomainError
        @test check_size(M, [1 2 3 4]) isa DomainError
        @test check_size(M, p, [1, 2]) isa DomainError
        @test check_size(M, p, [1 2 3 4]) isa DomainError
        @test embed(M, p) == p
        pE = similar(p)
        embed!(M, pE, p)
        @test pE == p
        P = [1.0 1.0 2.0]
        Q = similar(P)
        @test project!(M, Q, P) == project!(M, Q, P)
        @test project!(M, Q, P) == [1.0 1.0 0.0]

        XE = similar(X)
        embed!(M, XE, p, X)
        XE2 = embed(M, p, X)
        @test X == XE
        @test XE == XE2

        @test log(M, p, q) == q - p
        Y = similar(p)
        log!(M, Y, p, q)
        @test Y == q - p
        @test exp(M, p, X) == q
        r = similar(p)
        exp!(M, r, p, X)
        @test r == q
        @test distance(M, p, r) == norm(r - p)
    end

    @testset "AnotherPlaneManifold" begin
        M = AnotherPlaneManifold()
        p = [1.0, 2.0]
        pe = embed(M, p)
        @test pe == [1.0, 2.0, 0.0]
        X = [2.0, 3.0]
        Xe = embed(M, pe, X)
        @test Xe == [2.0, 3.0, 0.0]
        @test project(M, pe) == p
        @test project(M, pe, Xe) == X
    end

    @testset "Test nonimplemented fallbacks" begin
        @testset "Default Isometric Embedding Fallback Error Tests" begin
            M = NotImplementedEmbeddedManifold()
            A = zeros(2)
            # without any extra tests just the embedding is asked
            @test check_manifold_point(M, [1, 2]) === nothing
            @test check_tangent_vector(M, [1, 2], [3, 4]) === nothing
            @test norm(M, [1, 2], [2, 3]) ≈ sqrt(13)
            @test distance(M, [1, 2], [3, 4]) ≈ sqrt(8)
            @test inner(M, [1, 2], [2, 3], [2, 3]) ≈ 13
            @test_throws ErrorException manifold_dimension(M)
            # without any implementation the projections are the identity
            @test project(M, [1, 2]) == [1, 2]
            @test project(M, [1, 2], [2, 3]) == [2, 3]
            project!(M, A, [1, 2], [2, 3])
            @test A == [2, 3]
            @test vector_transport_direction(M, [1, 2], [2, 3], [3, 4]) == [2, 3]
            vector_transport_direction!(M, A, [1, 2], [2, 3], [3, 4])
            @test A == [2, 3]
            @test vector_transport_to(M, [1, 2], [2, 3], [3, 4]) == [2, 3]
            vector_transport_to!(M, A, [1, 2], [2, 3], [3, 4])
            @test A == [2, 3]
        end
        @testset "General Isometric Embedding Fallback Error Tests" begin
            M2 = NotImplementedEmbeddedManifold2()
            @test base_manifold(M2) == M2
            A = zeros(2)
            @test_throws ErrorException exp(M2, [1, 2], [2, 3])
            @test_throws ErrorException exp!(M2, A, [1, 2], [2, 3])
            @test_throws ErrorException log(M2, [1, 2], [2, 3])
            @test_throws ErrorException log!(M2, A, [1, 2], [2, 3])
            @test_throws ErrorException distance(M2, [1, 2], [2, 3])
            @test_throws ErrorException manifold_dimension(M2)
            @test_throws ErrorException project(M2, [1, 2])
            @test_throws ErrorException project!(M2, A, [1, 2])
            @test_throws ErrorException project(M2, [1, 2], [2, 3])
            @test_throws ErrorException project!(M2, A, [1, 2], [2, 3])
            @test_throws ErrorException vector_transport_along(M2, [1, 2], [2, 3], [[1, 2]])
            @test_throws ErrorException vector_transport_along(
                M2,
                [1, 2],
                [2, 3],
                [[1, 2]],
                ParallelTransport(),
            )
            @test_throws ErrorException vector_transport_along!(M2, A, [1, 2], [2, 3], [])
            @test_throws ErrorException vector_transport_direction(
                M2,
                [1, 2],
                [2, 3],
                [3, 4],
            )
            @test_throws ErrorException vector_transport_direction!(
                M2,
                A,
                [1, 2],
                [2, 3],
                [3, 4],
            )
            @test_throws ErrorException vector_transport_to(M2, [1, 2], [2, 3], [3, 4])
            @test_throws ErrorException vector_transport_to!(M2, A, [1, 2], [2, 3], [3, 4])
        end
        @testset "Nonisometric Embedding Fallback Error Rests" begin
            M3 = NotImplementedEmbeddedManifold3()
            @test_throws ErrorException inner(M3, [1, 2], [2, 3], [2, 3])
            @test_throws ErrorException manifold_dimension(M3)
            @test_throws ErrorException distance(M3, [1, 2], [2, 3])
            @test_throws ErrorException norm(M3, [1, 2], [2, 3])
            @test_throws ErrorException embed(M3, [1, 2], [2, 3])
            @test_throws ErrorException embed(M3, [1, 2])
        end
    end

    @testset "EmbeddedManifold decorator dispatch" begin
        TM = NotImplementedEmbeddedManifold() # transparently iso
        IM = NotImplementedEmbeddedManifold2() # iso
        AM = NotImplementedEmbeddedManifold3() # general
        for f in [
            embed,
            exp,
            get_basis,
            get_coordinates,
            get_vector,
            inverse_retract,
            log,
            norm,
            distance,
        ]
            @test ManifoldsBase.decorator_transparent_dispatch(f, AM) === Val(:parent)
        end
        for f in
            [project, retract, inverse_retract!, retract!, get_coordinates!, get_vector!]
            @test ManifoldsBase.decorator_transparent_dispatch(f, AM) === Val(:parent)
        end
        for f in [vector_transport_along, vector_transport_direction, vector_transport_to]
            @test ManifoldsBase.decorator_transparent_dispatch(f, AM) === Val(:parent)
        end
        for f in [mid_point, mid_point!]
            @test ManifoldsBase.decorator_transparent_dispatch(f, AM) === Val(:parent)
        end
        for f in [check_manifold_point, check_tangent_vector, exp!, inner, embed!]
            @test ManifoldsBase.decorator_transparent_dispatch(f, AM) ===
                  Val(:intransparent)
        end
        for f in [log!, manifold_dimension, project!]
            @test ManifoldsBase.decorator_transparent_dispatch(f, AM) ===
                  Val(:intransparent)
        end
        @test ManifoldsBase.decorator_transparent_dispatch(vector_transport_along!, AM) ===
              Val(:intransparent)
        @test ManifoldsBase.decorator_transparent_dispatch(
            vector_transport_to!,
            AM,
            1,
            1,
            1,
            1,
            1,
        ) === Val(:intransparent)
        @test ManifoldsBase.decorator_transparent_dispatch(
            vector_transport_direction!,
            AM,
        ) === Val(:parent)

        for f in [inner, norm]
            @test ManifoldsBase.decorator_transparent_dispatch(f, IM) === Val(:transparent)
        end
        for f in [inverse_retract!, retract!, mid_point!, distance]
            @test ManifoldsBase.decorator_transparent_dispatch(f, IM) === Val(:parent)
        end

        for f in [exp, inverse_retract, log, project, retract, mid_point, distance]
            @test ManifoldsBase.decorator_transparent_dispatch(f, TM) === Val(:transparent)
        end
        for f in [exp!, inverse_retract!, log!, project!, retract!, mid_point!]
            @test ManifoldsBase.decorator_transparent_dispatch(f, TM) === Val(:transparent)
        end
        for f in [vector_transport_along, vector_transport_direction, vector_transport_to]
            @test ManifoldsBase.decorator_transparent_dispatch(f, TM) === Val(:transparent)
        end
        for f in
            [vector_transport_along!, vector_transport_direction!, vector_transport_to!]
            @test ManifoldsBase.decorator_transparent_dispatch(f, TM) === Val(:transparent)
        end
        for t in [PoleLadderTransport(), SchildsLadderTransport()], M in [AM, TM, IM]
            @test ManifoldsBase.decorator_transparent_dispatch(
                vector_transport_to!,
                M,
                1,
                1,
                1,
                1,
                PoleLadderTransport(),
            ) === Val(:parent)
        end
        @test ManifoldsBase.decorator_transparent_dispatch(embed, TM) === Val{:parent}()
        @test ManifoldsBase.decorator_transparent_dispatch(embed!, TM) ===
              Val(:intransparent)
    end

    @testset "Explicit Embeddings using EmbeddedManifold" begin
        M = DefaultManifold(3, 3)
        N = DefaultManifold(4, 4)
        O = EmbeddedManifold(M, N)
        # first test with same length of sizes
        p = ones(3, 3)
        q = zeros(4, 4)
        qT = zeros(4, 4)
        qT[1:3, 1:3] .= 1.0
        embed!(O, q, p)
        @test norm(qT - q) == 0
        qM = embed(O, p)
        @test norm(project(O, qM) - p) == 0
        @test norm(qT - qM) == 0
        # test with different sizes, check that it only fills first element
        q2 = zeros(4, 4, 3)
        q2T = zeros(4, 4, 3)
        q2T[1:3, 1:3, 1] .= 1.0
        embed!(O, q2, p)
        @test norm(q2T - q2) == 0
        O2 = EmbeddedManifold(M, DefaultManifold(4, 4, 3))
        q2M = embed(O2, p)
        @test norm(q2T - q2M) == 0
        # wrong size error checks
        @test_throws DomainError embed!(O, zeros(3, 3), zeros(3, 3, 5))
        @test_throws DomainError embed!(O, zeros(3, 3), zeros(4, 4))
        @test_throws DomainError project!(O, zeros(3, 3, 5), zeros(3, 3))
        @test_throws DomainError project!(O, zeros(4, 4), zeros(3, 3))
        for f in [embed, project]
            @test ManifoldsBase.decorator_transparent_dispatch(f, O) === Val(:intransparent)
        end
    end
end


text:
## Types
#
    abstract type CompletionProvider end

    mutable struct NoCompletion <: CompletionProvider
    end

    mutable struct NormalCompletion <: CompletionProvider
        steps::Array{Function, 1}
        state::Int
        cmd::AbstractString
        itstart
        itend
        comp
        dotpos
        NormalCompletion() = new(Function[], 1, "", nothing, nothing, [""], -1:0)
    end

    mutable struct MethodCompletion <: CompletionProvider
        steps::Array{Function, 1}
        state::Int
        cmd::AbstractString
        itstart
        itend
        comp
        dotpos
        MethodCompletion() = new(Function[], 1, "", nothing, nothing, [""], -1:0)
    end

    mutable struct ArrayCompletion2 <: CompletionProvider
        steps::Array{Function, 1}
        state::Int
        cmd::AbstractString
        itstart
        itend
        comp
        dotpos
        indices
        ArrayCompletion2() = new(Function[], 1, "", nothing, nothing, [""], -1:0, Int[])
    end

    mutable struct TupleCompletion <: CompletionProvider
        steps::Array{Function, 1}
        state::Int
        cmd::AbstractString
        itstart
        itend
        comp
        func_names
        TupleCompletion() = new(Function[], 1, "", nothing, nothing, [""], [""])
    end

    mutable struct PathCompletion <: CompletionProvider
        steps::Array{Function, 1}
        state::Int
        cmd::AbstractString
        itstart
        itend
        comp
        dotpos
        PathCompletion() = new(Function[], 1, "", nothing, nothing, [""], -1:0)
    end

    mutable struct HistoryCompletion <: CompletionProvider
        steps::Array{Function, 1}
        state::Int
        cmd::AbstractString
        itstart
        itend
        comp
        dotpos
        HistoryCompletion() = new(Function[], 1, "", nothing, nothing, [""], -1:0)
    end

#end
#using CompletionProviders
##

function init_autocomplete(view::GtkTextView, t::EditorTab, replace=true; key=:tab)

    buffer = view.buffer[GtkTextBuffer]
    editor = parent(t)::Editor
    console = current_console(editor)

    #let's not autocomplete multiple lines
    (found, it_start, it_end) = selection_bounds(buffer)
    if found
        GtkTextUtils.nlines(it_start, it_end) > 1 && @goto exit
    end

    p = get_completion_provider(console, view, t, key)
    typeof(p) == NoCompletion && @goto exit

    if p.state <= length(p.steps)
        p.comp = p.steps[p.state]()
        p.state += 1
    else
        completions(p, t, completion_window.idx, console)
    end
    isempty(p.comp) && @goto exit

    if length(p.comp) == 1 && replace
        insert(p, formatcompletion(p, 1), buffer)
    else
        init_completion_window(view, p; mode = key)
    end

    return INTERRUPT

    @label exit
    visible(completion_window, false)
    return PROPAGATE
end

##

function get_completion_provider(console::Console, view::GtkTextView, t::EditorTab, key=:tab)
    buffer = view.buffer[GtkTextBuffer]
    it = get_text_iter_at_cursor(buffer)

    #change providers depending on which key we used
    pts = key == :tab ? [PathCompletion, ArrayCompletion2, NormalCompletion, MethodCompletion, TupleCompletion] : [HistoryCompletion]

    for pt in pts
        p = pt()
        select_text(p, console, buffer, it, t) && return p
    end
    return NoCompletion()
end

##

function test_completion_providers()
    t = get_current_tab()
    view = t.view
    p = get_completion_provider(console, view, t)

    typeof(p) == NoCompletion

    buffer = view.buffer[GtkTextBuffer]
    it = get_text_iter_at_cursor(buffer)

#    completions(p, t)
    init_autocomplete(view, t)
    p
end

#p = test_completion_providers()

## Default behavior

function formatcompletion(p::CompletionProvider, idx::Int)
    p.comp[idx]
end
function insert(p::CompletionProvider, s, buffer)
    replace_text(buffer, p.itstart, p.itend, s)
end

#####################
# Normal completion

function select_text(p::NormalCompletion, console, buffer, it, t)

    istextfile(t) && return false

    (cmd, its, ite) = select_word_backward(it, buffer, false)
    cmd = strip(cmd)
    isempty(cmd) && return false
    cmd[end] == '(' && return false #trying to complete a method

    cmd == "" && return false
    p.cmd = cmd
    p.itstart = its
    p.itend = ite
    true
end
function completions(p::NormalCompletion, t, idx, c::Console)
    if !isdefined(t, :autocomplete_words)
        t.autocomplete_words = [""]
    end
    comp, dotpos = extcompletions(p.cmd, t.autocomplete_words, c)
    p.comp = comp
    p.dotpos = dotpos
end

function formatcompletion(p::NormalCompletion, idx::Int)
    dotpos = p.dotpos.start
    prefix = dotpos > 1 ? p.cmd[1:dotpos-1] : ""
    prefix * p.comp[idx]
end

function insert(p::NormalCompletion, s, buffer)
    replace_text(buffer, p.itstart, p.itend, s)
end

#####################
# Methods completion

function select_text(p::MethodCompletion, console, buffer, it, t)
    istextfile(t) && return false

    (cmd, its, ite) = select_word_backward(it, buffer, false)
    cmd = strip(cmd)
    cmd == "" && return false
    cmd[end] != '(' && return false
    p.cmd = cmd
    p.itstart = its
    p.itend = ite
    true
end

function completions(p::MethodCompletion, t, idx, c::Console)
    # for some reason completion removes the first module when completing methods,
    # so I take note of it here and add it back at the end
    mods = split(p.cmd, '.')
    prefix = length(mods) > 2 ? mods[1] : ""

    comp, dotpos = GtkREPL.completions_in_module(p.cmd, c)
    dotpos = dotpos .+ (lastindex(prefix)-1)
    comp = [prefix != "" ? string(prefix, '.', c) : c for c in comp]

    p.comp = comp
    p.dotpos = dotpos
end

function formatcompletion(p::MethodCompletion, idx::Int)
    s = remove_filename_from_methods_def( p.comp[idx] )
    dotpos = p.dotpos.start
    prefix = dotpos > 1 ? p.cmd[1:dotpos-1] : ""
    prefix * s
end

function insert(p::MethodCompletion, s, buffer)
    replace_text(buffer, p.itstart, p.itend, s)
end

#####################
# Array completion

function select_text(p::ArrayCompletion2, console, buffer, it, t)

    istextfile(t) && return false

    (cmd, its, ite) = select_word_backward(it, buffer, false)
    cmd = strip(cmd)
    cmd == "" && return false
    length(cmd) == 1 && return false #we need something before [
    cmd[end] != '[' && return false

    p.cmd = cmd[1:end-1]
    p.itstart = its
    p.itend = ite
    true
end

function compact_output(x)
    io = IOBuffer()
    show(IOContext(io, :compact => true, :limit => true, :displaysize => (8,20)), "text/plain", x)
    String(take!(io))
end

function completions(p::ArrayCompletion2, t, idx, c::Console)

    var = Symbol(p.cmd)
    #check if variable is defined
    !remotecall_fetch(isdefined, worker(c), c.eval_in, var) && return false

    ex =  :( (size = size($var), type = typeof($var)) )
    s = remotecall_fetch(Core.eval, worker(c), c.eval_in, ex)

    elems = [remotecall_fetch(Core.eval, worker(c), c.eval_in, Expr(:ref, var, i)) for i = 1:min(sum(s.size),20)]
    elems = compact_output.(elems)
    elems = [replace(e, '\n' => ' ') for e in elems]

    idx   = [CartesianIndices(s.size)[i].I for i = 1:min(sum(s.size),20)]
    elems = [string(idx, ": ", el) for (idx,el) in zip(idx,elems)]

    comp  = ["$(s.size) - $(s.type)"; string.(elems)]
    #comp = (s, string.(elems))
    p.indices = idx
    p.comp = comp
end

#TODO insert index into the brackets
function insert(p::ArrayCompletion2, s, buffer)
    replace_text(buffer, p.itend, p.itend, s)
end

function formatcompletion(p::ArrayCompletion2, idx::Int)
    idx == 1 && return "" # the first element of comp shows the size of the array
    i = p.indices[idx-1]
    string( join(i, ','), ']')
end

#####################
#Tuple completion

function select_text(p::TupleCompletion, console, buffer, it, t)
    istextfile(t) && return false

    (found, tu, itstart) = select_tuple(it, buffer)
    !found && return false

    p.cmd = tu
    p.itstart = itstart
    p.itend = itstart
    true
end

function completions(p::TupleCompletion, t, idx, c::Console)

    args = tuple_to_types(p.cmd, c)
    isempty(args) && return
    m = methods_with_tuple(args)
    comp = map(string, m)
    func_names = [string(x.name) for x in m]
    p.comp = comp
    p.func_names = func_names
end

function formatcompletion(p::TupleCompletion, idx::Int)
    p.func_names[idx]
end

function insert(p::TupleCompletion, s, buffer)
    insert!(buffer, p.itstart, s)
end

#####################
# PathCompletion

function formatcompletion(p::PathCompletion, idx::Int)
    dotpos = p.dotpos.start
    cmd = p.cmd[2:end]#remove the leading "
    prefix = dotpos > 1 ? cmd[1:dotpos-1] : ""

    '"' * prefix * p.comp[idx]
end

function select_text(p::PathCompletion, console, buffer, it, t)
    istextfile(t) && return false

    (cmd, its, ite) = get_current_line_text(buffer)
    ite = mutable(get_text_iter_at_cursor(buffer))
    cmd = (its:ite).text[String]

    cmd == "" && return false
    idx = findlast(c->c=='"', cmd)

    idx == nothing && return false
    idx == length(cmd) && return false
    its += idx-1
    cmd = cmd[idx:end]

    #check if it's a folder
    cmd_s = strip(cmd[2:end])
    comp, dotpos = remotecall_fetch(REPL.REPLCompletions.complete_path, GtkIDE.worker(console), cmd_s, lastindex(cmd_s))
    #if not continue with normal completion
    isempty(comp) && return false

    p.cmd = cmd
    p.itstart = its
    p.itend = ite
    true
end

function completions(p::PathCompletion, t, idx, c::Console)
    cmd = p.cmd
    cmd = strip(cmd[2:end])#remove the leading "
    comp, dotpos = remotecall_fetch(REPL.REPLCompletions.complete_path, worker(c), cmd, lastindex(cmd))
    comp = [c.path for c in comp]

    p.dotpos = dotpos
    p.comp = comp
end

#####################
# HistoryCompletion

function select_text(p::HistoryCompletion, console, buffer, it, t)

    istextfile(t) && return false

    (found, its, ite) = selection_bounds(buffer)
    if found
        cmd = (its:ite).text[String]
    else
        #(cmd, its, ite) = select_word_backward(it, buffer, false)
        (cmd, its, ite) = get_current_line_text(buffer)
    end
    cmd = strip(cmd)
    isempty(cmd) && return false

    cmd == "" && return false
    p.cmd = cmd
    p.itstart = its
    p.itend = ite
    true
end

function completions(p::HistoryCompletion, t, idx, c::Console)

    cmd = strip(p.cmd)
    idx = GtkREPL.search(c.history, cmd)

    p.dotpos = -1
    p.comp = unique(c.history.history[idx])
end

function select_text(p::HistoryCompletion, console, buffer, it, t)

    istextfile(t) && return false

    (found, its, ite) = selection_bounds(buffer)
    if found
        cmd = (its:ite).text[String]
    else
        #(cmd, its, ite) = select_word_backward(it, buffer, false)
        # select the whole line, so it works like in the console
        (cmd, its, ite) = get_current_line_text(buffer)
    end
    cmd = strip(cmd)
    isempty(cmd) && return false

    cmd == "" && return false
    p.cmd = cmd
    p.itstart = its
    p.itend = ite
    true
end


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Heerden2013")

import CSV
import MAT
using DataFrames
using Serialization

## -------------------------------------------------------------------
# Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in the Julia Pkg REPL to install the
# package, then you must activate the package enviroment (see README)
import Chemostat_Heerden2013
const Hd  = Chemostat_Heerden2013.HeerdenData;
const ECC = Chemostat_Heerden2013.EColiCore

## -------------------------------------------------------------------
# run add "https://github.com/josePereiro/Chemostat" in the
# julia Pkg REPL for installing the package
import Chemostat
const ChU = Chemostat.Utils
const ChSS = Chemostat.SteadyState
const ChLP = Chemostat.LP

## -------------------------------------------------------------------
# Tools
function partial_test(model, title  = "PARTIAL TEST")
    fbaout = ChLP.fba(model, ECC.BIOMASS_IDER);
    ChU.tagprintln_inmw(title,
        "\nsize:             ", size(model),
        "\nobj_ider:         " , ECC.BIOMASS_IDER,
        "\nfba obj_val:      ", ChU.av(model, fbaout, ECC.BIOMASS_IDER),
        "\nmax exp obj_val:  ", maximum(Hd.val("D")),
    )
end

## -------------------------------------------------------------------
# BASE MODEL
# Load Mat file
println("Original .mat model")
src_file = ECC.MODEL_RAW_MAT_FILE
mat_model = MAT.matread(src_file)["e_coli_core"]
model = ChU.MetNet(mat_model; reshape = true);

## -------------------------------------------------------------------
# the maximal experimental growth rate in Heerden2013 is ~0.2 1/h
# The raw model present a growth rate bigger than that, so it is ok
# to use it directly as base model
partial_test(model)

## -------------------------------------------------------------------
# Set bounds
# The abs maximum bounds will be set to ECC.MAX_ABS_BOUND
ChU.tagprintln_inmw("CLAMP BOUNDS",
    "\nabs max bound: ", ECC.MAX_ABS_BOUND
)
ChU.clampfields!(model, [:lb, :ub];
    abs_max = ECC.MAX_ABS_BOUND, zeroth = 0.0)
partial_test(model)

## -------------------------------------------------------------------
# Exchanges
exchs = ChU.exchanges(model)
ChU.tagprintln_inmw("SETTING EXCHANGES")
# To control the intakes just the metabolites defined in the
# base_intake_info (The minimum medium) will be opened.
# The base model will be constraint as in a cultivation with xi = 1.0
# see Cossios paper (see README)
ξ = 1.0
# println("Minimum medium: ", ECC.base_intake_info)
foreach(exchs) do idx
    ChU.ub!(model, idx, ECC.MAX_ABS_BOUND) # Opening all outakes
    ChU.lb!(model, idx, 0.0) # Closing all intakes
end

# see Cossios paper (see README) for details in the Chemostat bound constraint
ChSS.apply_bound!(model, ξ, ECC.base_intake_info);

## -------------------------------------------------------------------
# Eliminate similars
ChU.tagprintln_inmw("ELIMINATE SIMILAR RXNS")
ChU.summary.([model], ["SUCDi", "FRD7"])
ChU.bounds!(model, "SUCDi", 0.0, 0.0)
ChU.bounds!(model, "FRD7", -ECC.MAX_ABS_BOUND, ECC.MAX_ABS_BOUND)
ChU.summary.([model], ["SUCDi", "FRD7"])

## -------------------------------------------------------------------
# Exch_met_map

# A quick way to get exchages from mets and te other way around
exch_met_map = Dict()
for exch_ in model.rxns[findall((id) -> any(startswith.(id, ["EX_", "DM_"])), model.rxns)]
    mets_ = model.mets[ChU.rxn_mets(model, exch_)]
    length(mets_) != 1 && continue
    exch_met_map[exch_] = mets_[1]
    exch_met_map[mets_[1]] = exch_
end;

# saving
ChU.save_data(ECC.EXCH_MET_MAP_FILE, exch_met_map)

## -------------------------------------------------------------------
partial_test(model)

# FVA PREPROCESSING
model = ChLP.fva_preprocess(model,
    batchlen = 1,
    check_obj = ECC.BIOMASS_IDER,
    verbose = true
)
partial_test(model)

## -------------------------------------------------------------------
ChU.summary(model)

## -------------------------------------------------------------------
# Saving model
dict_model = model |> ChU.struct_to_dict
ChU.save_data(ECC.BASE_MODELS_FILE, dict_model);


text:
"""
A "type" to indicate a state of an automaton
"""
State = String

"""
An automaton having:
 - A set of start states q₀
 - A set of final/accepting states F
 - A set of input symbols (alphabet) Σ
 - A set of states Q
 - A transition function δ

Notations are from Wikipedia:
https://en.wikipedia.org/wiki/Deterministic_finite_automaton
"""
struct Automaton
    q₀::State
    F::Set{State}
    Σ::Set{Char}
    Q::Set{State}
    δ::Dict{Tuple{State,Char},Vector{State}}
end

function Automaton(δ::AbstractDict{},
                   q₀::State,
                   F::AbstractSet{State})
    # infer alphabet from transitions
    Σ = Set{Char}()

    # Q will be set of all states
    Q = Set{State}()
    # reached will be set of destination states
    reached = Set{State}()
    for (k, v) in δ
        # destination states
        for d_state in v
            push!(Q, State(d_state))
        end

        # non-final source state
        push!(Q, State(k[1]))

        push!(Σ, k[2])  # add to alphabet
    end

    Automaton(q₀, F, Σ, Q, δ)
end

abstract type AbstractFiniteAutomaton end

"""
    evaluate!(M::AbstractFiniteAutomaton, w::String)

Evaluates the input string using the given automaton.
The automaton is reset before evaluating.

Returns `true` if string is accepted and
`false` otherwise.
"""
function evaluate!(M::AbstractFiniteAutomaton, w)
    reset!(M)

    for a in w
        transition!(M, a) || return false
    end

    accepting(M)
end

max_issues_repo_name: MadhavJivrajani/DFA-Sim

max_forks_repo_name: MadhavJivrajani/DFA-Sim


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

include("../helpers/TikzUtils.jl")
include("../helpers/save_helpers.jl")

using RunCases
using Obj2Dict

using TikzPictures
import PGFPlots: Plots,Axis

function irs_2ac_mcts_script(encounter::Int64, iterations::Vector{Int64})
  #for 1 encounter

  cases = generate_cases(("sim_params.encounter_number",
                          [encounter]),
                         ("mcts_params.maxtime_s",
                          [realmax(Float64)]),
                         ("mcts_params.n",
                          iterations))

  add_field!(cases, "study_params.fileroot", n -> "trajSaveMCTS_enc$(encounter)_n$(n)",
             ["mcts_params.n"])

  trajSave(MCTSStudy(), cases)

end

function irs_vis{T<:AbstractString}(files::Vector{T}; outfileroot::AbstractString="irs_vis")

  #attributes defined by callbacks
  #[encounter number, iterations, total reward, isnmac]
  callbacks = [x -> sv_encounter_id(x)[1], sv_dpw_iterations, sv_reward, sv_nmac]

  #encounter files (rows) x attributes (cols)
  data = Array(Any, length(files), length(callbacks))

  for (row, file) in enumerate(files)
    d = trajLoad(file)
    data[row, :] = map(f -> f(d), callbacks)
  end

  td = TikzDocument()
  plotArray = Plots.Plot[]

  for enc in unique(data[:, 1])
    rows = find(x -> x[1] == enc, data[:, 1])
    sort!(rows, by=r -> data[r, 2])
    xs = convert(Vector{Float64}, data[rows, 2])
    ys = convert(Vector{Float64}, data[rows, 3])

    push!(plotArray,Plots.Linear(xs, ys, legendentry = "Encounter $(enc)"))
  end

  tp = PGFPlots.plot(Axis(plotArray,
                          xlabel = "MCTS Iterations",
                          ylabel = "Total Reward",
                          title = "MCTS DPW Iterations vs. Total Reward",
                          legendPos = "south east"))

  cap = "MCTS DPW Iterations vs. Total Reward"

  add_to_document!(td, tp, cap)

  outfile = string(outfileroot, ".pdf")
  TikzPictures.save(PDF(outfile), td)
  outfile = string(outfileroot, ".tex")
  TikzPictures.save(TEX(outfile), td)

end

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
module MaximumGeneratorProfitTCRDD

#  Packages
using Cbc
using D3TypeTrees
using DataFrames
using Dates
using DataStructures
using GLPK
using Ipopt
using JuMP
using LinearAlgebra
using Plots
using PowerSimulations
using PowerSystems
using TimeSeries

#  Package Constants
const PSY = PowerSystems
const PSI = PowerSimulations

#  Functions
# Power System Auxiliary Functions for TCRDD Algorithm
include("powersystemaux/assign_lines_index.jl")
include("powersystemaux/change_slack!.jl")
include("powersystemaux/get_binding_lines.jl")
#include("powersystemaux/get_Pg_thermal_bindingPg.jl")
#include("powersystemaux/get_Pg_thermal_constprice.jl")
#include("powersystemaux/get_Pg_thermal_nonzeroslope.jl")
#include("powersystemaux/get_Pg_thermal_slack.jl")
include("powersystemaux/get_PTDF_bindingLines.jl")
include("powersystemaux/get_PTDF_load.jl")
include("powersystemaux/get_PTDF_thermal_bindingPg.jl")
include("powersystemaux/get_PTDF_thermal_constprice.jl")
include("powersystemaux/get_PTDF_thermal_nonzeroslope.jl")
include("powersystemaux/get_PTDF_thermal_slack.jl")
#include("powersystemaux/get_thermal_bindingPg.jl")
#include("powersystemaux/get_thermal_constprice.jl")
#include("powersystemaux/get_thermal_nonzeroslope.jl")
include("powersystemaux/get_thermal_slack.jl")
include("powersystemaux/opf_PTDF.jl")

# Functions for TCRDD Algorithm

include("algorithm/bid_argmax_profit.jl")
include("algorithm/bisection_loop!.jl")
include("algorithm/compare_tan_inv_RDC.jl")
include("algorithm/d2f_PgPg_nonzeroslope.jl")
include("algorithm/disp_plots.jl")
include("algorithm/f_TCRDD.jl")
include("algorithm/residual_demand.jl")

# Main Functions
include("maxGenProfit_tcrdd.jl")
include("plots_TCRDD.jl")

# Type

# Accessor functions

# Useful Calculation Functions
#export change_slack!

# Main Functions
export maxGenProfit_tcrdd
export plots_TCRDD

end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
include("../src/RationalFunctionFields.jl")

using .RationalFunctionFields
using .RationalFunctionFields: load_generators, groebner_ideal,
                            RationalFunctionField, compute_groebner!, gens

Singular = RationalFunctionFields.Singular
Logging = RationalFunctionFields.Logging

logger = Logging.SimpleLogger(stderr, Logging.Debug)
Logging.global_logger(logger)

set = load_generators("../data/Goodwin.txt")

println(set, "\n", length(set))

FF = RationalFunctionField(set)

compute_groebner!(FF)

println( "\nnew gens ", gens(FF) )


text:
function palign(
    lseq,
    J::Array{Float64,4},
    h::Array{Float64,2},
    λo::Array{Float64,1},
    λe::Array{Float64,1};
    verbose::Bool = true,
    maxiter::Integer = 100,
    nprint::Integer = 500,
    epsconv::Real = 1e-4,
    seed::Int = 0,
    damp::Real = 0,
    μint::Real = 1.0,
    μext::Real = 1.0,
    λ::Real = 1.0,
    T0::Bool = false,
    mindec::Int = 10,
    ctype::Symbol = "amino",
)

    Random.seed!(seed)
    en = Inf
    seq = Seq(lseq.header, lseq.strseq, lseq.intseq, ctype) # crappy way to avoid recompilation
    jh = Jh(J, h)
    alg =
        Alg(verbose, maxiter, epsconv, mindec, damp, μext, λo, λe, μint, nprint)
    pbf = PBF(jh, seq)
    aux = lseq.strseq
    allvar = AllVar(pbf, jh, seq, alg)
    @extract pbf:L N
    @extract pbf:newP
    if T0 == false
        if verbose
            println("Run Belief Propagation, large-connectivity approximation, for proteins alignment")
            println("L $L = length of the Potts model")
            println("N $N = length of the sequence A to be aligned")
            println("A: $aux")
        end
        maxiter, flagconv, en = update!(allvar)
    else
        if verbose
            println("Run Belief Propagation, large-connectivity approximation, for proteins alignment")
            println("L $L = length of the Potts model")
            println("N $N = length of the sequence A to be aligned")
            println("A: $aux")
            println("T = 0, using β → +∞ limit")
        end
        maxiter, flagconv, en = updateT0!(allvar)
    end
    return maxiter, flagconv, allvar, en
end

function palign(
    filefasta::String,
    J::Array{Float64,4},
    h::Array{Float64,2},
    λo::Array{Float64,1},
    λe::Array{Float64,1},
    ctype::Symbol;
    kwdargs...,
)
    seq = readunalignedfasta(filefasta, ctype)[1]
    palign(seq, J, h; kwdargs...)
end


text:
@testset "@cunion" begin
	@cunion OpaqueUnion
	@cunion OpaqueUnion
	@test_throws ErrorException sizeof(OpaqueUnion)
	@test sizeof(Ptr{OpaqueUnion}) == sizeof(Ptr{Cvoid})

	@cunion EmptyUnion {
	} __packed__
	@test sizeof(EmptyUnion) == 0

	@cunion CcharUnion {
		c::Cchar
	} __packed__
	@test sizeof(CcharUnion) == sizeof(Cchar)
	@test :c in propertynames(CcharUnion)

	@cunion CintUnion {
		i::Cint
	} __packed__
	@test sizeof(CintUnion) == sizeof(Cint)
	@test :i in propertynames(CintUnion)

	@cunion CunionUnion {
		i::CintUnion
	} __packed__
	@test sizeof(CunionUnion) == sizeof(CintUnion)
	@test :i in propertynames(CunionUnion)

	@cunion CintCcharUnion {
		i::CintUnion
		c::CcharUnion
	} __packed__
	@test sizeof(CintCcharUnion) == max(sizeof(CintUnion), sizeof(CcharUnion))
	@test :i in propertynames(CintCcharUnion)
	@test :c in propertynames(CintCcharUnion)

	@cunion CcharCintPackedUnion {
		c::CcharUnion
		i::CintUnion
	} __packed__
	@test sizeof(CcharCintPackedUnion) == max(sizeof(CintUnion), sizeof(CcharUnion))
	@test :c in propertynames(CcharCintPackedUnion)
	@test :i in propertynames(CcharCintPackedUnion)

	@cunion CcharCintAlignedUnion {
		c::CcharUnion
		@calign sizeof(CintUnion)
		i::CintUnion
	} __packed__
	@test sizeof(CcharCintAlignedUnion) == max(max(sizeof(CintUnion), sizeof(Cint)), max(sizeof(CcharUnion), sizeof(Cint)))
	@test :c in propertynames(CcharCintAlignedUnion)
	@test :i in propertynames(CcharCintAlignedUnion)

	@cunion CunionUnion {
		@cunion {
			c::Cuchar[4]
			s::Cushort[2]
			i::Cuint
		} __packed__
	} __packed__
	@test :c in propertynames(CunionUnion)
	@test :s in propertynames(CunionUnion)
	@test :i in propertynames(CunionUnion)

	@cunion PtrToUnionUnion
	@cunion PtrToPtrToUnionUnion {
		p::Ptr{PtrToUnionUnion}
	} __packed__
	@cunion PtrToUnionUnion {
		p::Ptr{PtrToUnionUnion}
	} __packed__
	@test sizeof(PtrToPtrToUnionUnion) == sizeof(Ptr)
	@test :p in propertynames(PtrToPtrToUnionUnion)
	@test sizeof(PtrToUnionUnion) == sizeof(Ptr)
	@test :p in propertynames(PtrToUnionUnion)

	cuu = CunionUnion()
	@test cuu.i == 0
	@test cuu.c[1] == 0
	@test cuu.c[4] == 0
	@test cuu.s[1] == 0
	@test cuu.s[2] == 0
	cuu.i = 0xff0000ff
	@test cuu.i == 0xff0000ff
	@test cuu.c[1] == 0xff
	@test cuu.c[4] == 0xff
	@test cuu.s[1] != 0
	@test cuu.s[2] != 0

	@cunion Cint32BitfieldUnion {
		(i:32)::Cint
	} __packed__
	@test sizeof(Cint32BitfieldUnion) == sizeof(Cint)
	@test :i in propertynames(Cint32BitfieldUnion)

	@cunion Cuint32BitfieldUnion {
		(u:32)::Cuint
	} __packed__
	@test sizeof(Cuint32BitfieldUnion) == sizeof(Cuint)
	@test :u in propertynames(Cuint32BitfieldUnion)

	@cunion Cint2BitfieldUnion {
		(i:2)::Cint
	} __packed__
	@test sizeof(Cint2BitfieldUnion) == 1
	@test :i in propertynames(Cint2BitfieldUnion)

	@cunion Cuint2BitfieldUnion {
		(u:2)::Cuint
	} __packed__
	@test sizeof(Cuint2BitfieldUnion) == 1
	@test :u in propertynames(Cuint2BitfieldUnion)

	@cunion Cuint32Cint32BitfieldsUnion {
		(u:32)::Cuint
		(i:32)::Cint
	} __packed__
	@test sizeof(Cuint32Cint32BitfieldsUnion) == sizeof(Cuint)
	@test :u in propertynames(Cuint32Cint32BitfieldsUnion)
	@test :i in propertynames(Cuint32Cint32BitfieldsUnion)

	@cunion Cuint16Cint16BitfieldsUnion {
		(u:16)::Cuint
		(i:16)::Cint
	} __packed__
	@test sizeof(Cuint16Cint16BitfieldsUnion) == 2
	@test :u in propertynames(Cuint16Cint16BitfieldsUnion)
	@test :i in propertynames(Cuint16Cint16BitfieldsUnion)

	@cunion Cint16Cint16BitfieldsUnion {
		(i1:16)::Cint
		(i2:16)::Cint
	} __packed__
	@test sizeof(Cint16Cint16BitfieldsUnion) == 2

	@cunion Cint16CalignCint16BitfieldsUnion {
		(i1:16)::Cint
		@calign sizeof(Cint)
		(i2:16)::Cint
	} __packed__
	@test sizeof(Cint16CalignCint16BitfieldsUnion) == sizeof(Cint)

	@cunion Cint16CintBitfieldsUnion {
		(bf:16)::Cint
		i::Cint
	} __packed__
	@test sizeof(Cint16CintBitfieldsUnion) == sizeof(Cint)

	bfUnion = Cuint16Cint16BitfieldsUnion()
	@test bfUnion.u == 0
	@test bfUnion.i == 0
	bfUnion.i = Cint(-1)
	@test bfUnion.u == 0xffff
	@test bfUnion.i == -1
	bfUnion.u = Cuint(0x0001)
	@test bfUnion.u == 0x0001
	@test bfUnion.i == 1

	bfUnion = Cint16CalignCint16BitfieldsUnion()
	@test bfUnion.i1 == 0
	@test bfUnion.i2 == 0
	bfUnion.i1 = -1
	@test bfUnion.i1 == -1
	@test bfUnion.i2 == -1
	bfUnion.i2 = 0
	@test bfUnion.i1 == 0
	@test bfUnion.i2 == 0
	bfUnion.i2 = 32768
	@test bfUnion.i1 < 0
	@test bfUnion.i2 < 0
end


text:
using SettingsBundles
using Base.Test

# write your own tests here
@test 1 == 1

settings = SettingsBundle()
add!(settings, SettingsYAMLFile("/Users/dustinrb/.julia/v0.4/SettingsBundles/test/test.yaml"))
add!(settings, SettingsYAMLFile("/Users/dustinrb/.julia/v0.4/SettingsBundles/test/test.json"))
settings["something"]


text:
# CUDNN_RNN_RELU: Stock RNN with ReLu activation
# CUDNN_RNN_TANH: Stock RNN with tanh activation
# CUDNN_LSTM:     LSTM with no peephole connections
# CUDNN_GRU:      Using h' = tanh(r * Uh(t-1) + Wx) and h = (1 - z) * h' + z * h(t-1)

# param layout:
# RNN: [weight, bias] × [input, hidden]
# GRU: [weight, bias] × [input, hidden] × [reset, update, newmem]
# LSTM: [weight, bias] × [input, hidden] × [input, forget, newmem, output]

using LinearAlgebra

function params(w::CuVector, input, hidden, n = 1)
  slice(offset, shape) = reshape(view(w, offset.+(1:prod(shape))), shape)
  wx = slice(0, (input, hidden*n))
  wh = slice(length(wx), (hidden, hidden*n))
  bias = view(w, length(wx)+length(wh) .+ (1:hidden*n))
  (wx, wh), bias
end

mutable struct RNNDesc{T}
  mode::cudnnRNNMode_t
  input::Int
  hidden::Int
  params::CuVector{T}
  weights::NTuple{2,CuMatrix{T}}
  bias::CuVector{T}
  ptr::Ptr{Nothing}
end

Base.unsafe_convert(::Type{Ptr{Nothing}}, d::RNNDesc) = d.ptr

function rnnParamSize(T, r, input)
  size = Csize_t[0]
  cudnnGetRNNParamsSize(handle(), r, TensorDesc(T, (1,input,1)), size, cudnnDataType(T))
  return Int(size[])÷sizeof(T)
end

ngates(mode) = [1, 1, 4, 3][mode+1]
ngates(r::RNNDesc) = ngates(r.mode)

function RNNDesc{T}(mode::cudnnRNNMode_t, input::Int, hidden::Int; layers = 1) where T
  d = [C_NULL]
  cudnnCreateRNNDescriptor(d)

  dropoutDesc = DropoutDesc(0)
  inputMode = CUDNN_LINEAR_INPUT
  direction = CUDNN_UNIDIRECTIONAL
  algo = CUDNN_RNN_ALGO_STANDARD
  cudnnSetRNNDescriptor_v6(handle(),d[],hidden,layers,dropoutDesc,inputMode,direction,mode,algo,cudnnDataType(T))

  w = CUDA.zeros(T, rnnParamSize(T, d[], input))
  # TODO: avoid reserve allocation here
  rd = RNNDesc{T}(mode, input, hidden, w, params(w, input, hidden, ngates(mode))..., d[])
  finalizer(rd) do x
    cudnnDestroyRNNDescriptor(x)
  end
  return rd
end

function setweights!(d::RNNDesc, Wi, Wh, b)
  transpose!(d.weights[1], Wi)
  transpose!(d.weights[2], Wh)
  copyto!(d.bias, b)
  return
end

function cudnnGetRNNTrainingReserveSize(r::RNNDesc, seqlen, xdesc)
  size = Csize_t[0]
  cudnnGetRNNTrainingReserveSize(handle(), r, seqlen, xdesc, size)
  return Int(size[])
end

function cudnnRNNForward(rnn::RNNDesc{T}, seqlen, xd, x, hd, h, cd, c, wd, w, yd, y, hod,
                         ho, cod, co, reserve=nothing) where T
  @workspace size=@argout(
      cudnnGetRNNWorkspaceSize(handle(), rnn, seqlen, xd,
                               out(Ref{Csize_t}()))
    )[] workspace->begin
      if reserve == nothing
        cudnnRNNForwardInference(handle(), rnn, seqlen, xd, x, hd, h, cd, c, wd, w, yd, y,
                                hod, ho, cod, co, workspace, sizeof(workspace))
      else
        cudnnRNNForwardTraining(handle(), rnn, seqlen, xd, x, hd, h, cd, c, wd, w, yd, y,
                                hod, ho, cod, co, workspace, sizeof(workspace),
                                reserve, sizeof(reserve))
      end
    end
end

xDesc(x) = [TensorDesc(eltype(x), (1, size(x, 1), size(x, 2)))]

hDesc(h::Nothing) = C_NULL, CU_NULL
hDesc(x::Integer) = (@assert x == 0; hDesc(nothing))
function hDesc(h::CuArray)
  TensorDesc(eltype(h), (size(h, 1), size(h, 2), 1)), h
end

# TODO: can we just manipulate strides here?
# TODO: should use repmat, but this isn't implemented.
hBatch(x::AbstractVector, h::CuVector) = h
hBatch(x::AbstractMatrix, h::CuVector) = h .* CUDA.ones(1, size(x, 2))
hBatch(x::AbstractMatrix, h::CuMatrix) = h .* CUDA.ones(1, size(h,2) == 1 ? size(x,2) : 1)

function forward(rnn::RNNDesc{T}, x::CuArray{T}, h_::CuArray{T}, c_ = nothing, train = Val{false}) where T
  h = hBatch(x, h_)
  c = c_ == nothing ? nothing : hBatch(x, c_)
  @assert size(x, 1) == rnn.input
  @assert size(h, 1) == rnn.hidden
  @assert size(x, 2) == size(h, 2)
  seqLength = 1
  xdesc = xDesc(x)
  y = x isa AbstractVector ? similar(x, rnn.hidden) : similar(x, rnn.hidden, size(x, 2))
  ho = similar(h)
  ydesc = xDesc(y)
  reserve = train == Val{true} ?
    CuVector{UInt8}(undef, cudnnGetRNNTrainingReserveSize(rnn, seqLength, xdesc)) :
    nothing
  co = c == nothing ? c : similar(c)
  cudnnRNNForward(rnn, seqLength,
                  xdesc, x,
                  hDesc(h)...,
                  hDesc(c)...,
                  FilterDesc(T, (1, 1, length(rnn.params))), rnn.params,
                  ydesc, y,
                  hDesc(ho)...,
                  hDesc(co)...,
                  reserve)
  result = c == nothing ? (y, ho) : (y, ho, co)
  return train == Val{true} ? (reserve, result) : result
end

forwardTrain(rnn::RNNDesc{T}, x::CuArray{T}, h::CuArray{T}, c = nothing) where T =
  forward(rnn, x, h, c, Val{true})

function cudnnRNNBackwardData(rnnDesc, seqLength, yDesc, y, dyDesc, dy, dhyDesc,
                              dhy, dcyDesc, dcy, wDesc, w, hxDesc, hx, cxDesc, cx, dxDesc,
                              dx, dhxDesc, dhx, dcxDesc, dcx, reserve)
  @workspace size=@argout(
      cudnnGetRNNWorkspaceSize(handle(), rnnDesc, seqLength, dxDesc,
                               out(Ref{Csize_t}()))
    )[] workspace->begin
      cudnnRNNBackwardData(handle(), rnnDesc, seqLength, yDesc, y, dyDesc, dy, dhyDesc,
                           dhy, dcyDesc, dcy, wDesc, w, hxDesc, hx, cxDesc, cx, dxDesc,
                           dx, dhxDesc, dhx, dcxDesc, dcx, workspace, sizeof(workspace),
                           reserve, sizeof(reserve))
    end
end

function backwardData(rnn::RNNDesc{T}, y, dy_, dho, dco, h, c, reserve) where T
  # Same as above, any more efficient way?
  dy = dy_ isa Integer ? zero(y) : dy_
  yd = xDesc(y)
  dx = y isa AbstractVector ? similar(dy, rnn.input) : similar(dy, rnn.input, size(dy, 2))
  dh = similar(h)
  dc = c == nothing ? nothing : similar(c)
  cudnnRNNBackwardData(rnn, 1, yd, y, yd, dy, hDesc(dho)..., hDesc(dco)...,
                       FilterDesc(T, (1, 1, length(rnn.params))), rnn.params, hDesc(h)...,
                       hDesc(c)..., xDesc(dx), dx, hDesc(dh)..., hDesc(dc)..., reserve)
  return c == nothing ? (dx, dh) : (dx, dh, dc)
end

backwardData(rnn, y, dy, dho, hx, reserve) =
  backwardData(rnn, y, dy, dho, nothing, hx, nothing, reserve)

function cudnnRNNBackwardWeights(rnnDesc, seqLength, xDesc, x, hxDesc, hx, yDesc,
                                 y, dwDesc, dw, reserve)
  @workspace size=@argout(
      cudnnGetRNNWorkspaceSize(handle(), rnnDesc, seqLength, xDesc,
                               out(Ref{Csize_t}()))
    )[] workspace->begin
      cudnnRNNBackwardWeights(handle(), rnnDesc, seqLength, xDesc, x, hxDesc, hx, yDesc,
                              y, workspace, sizeof(workspace), dwDesc, dw,
                              reserve, sizeof(reserve))
    end
end

function backwardWeights(rnn::RNNDesc{T}, x, h, y, reserve) where T
  dw = zero(rnn.params)
  cudnnRNNBackwardWeights(rnn, 1, xDesc(x), x, hDesc(h)..., xDesc(y), y,
                          FilterDesc(T, (1, 1, length(dw))), dw, reserve)
  return params(dw, rnn.input, rnn.hidden, ngates(rnn))
end

function pullback(rnn::RNNDesc{T}, x::CuArray{T}, h::CuArray{T}) where T <: Union{Float32,Float64}
  reserve, (y, ho) = CUDNN.forwardTrain(rnn, x, h)
  return (y, ho), function (dy, dho)
    h_ = CUDNN.hBatch(x, h)
    dx, dh = CUDNN.backwardData(rnn, y, dy, dho, h_, reserve)
    (dWi, dWh), db = CUDNN.backwardWeights(rnn, x, h_, y, reserve)
    return (x = dx, h = dh, Wi = dWi, Wh = dWh, b = db)
  end
end

function pullback(rnn::RNNDesc{T}, x::CuArray{T}, h::CuArray{T}, c::CuArray{T}) where T <: Union{Float32,Float64}
  reserve, (y, ho, co) = CUDNN.forwardTrain(rnn, x, h, c)
  return (y, ho, co), function (dy, dho, dco)
    h_ = CUDNN.hBatch(x, h)
    c_ = CUDNN.hBatch(x, c)
    dx, dh, dc = CUDNN.backwardData(rnn, y, dy, dho, dco, h_, c_, reserve)
    (dWi, dWh), db = CUDNN.backwardWeights(rnn, x, h_, y, reserve)
    return (x = dx, h = dh, c = dc, Wi = dWi, Wh = dWh, b = db)
  end
end


text:
module TimeSync

using Dates

module Lib
    using CxxWrap, Libdl

    const DEPSDIR = joinpath(dirname(@__DIR__), "deps")
    Libdl.dlopen(joinpath(DEPSDIR, "libsync.so"))

    @wrapmodule(joinpath(DEPSDIR, "libsync.so"))

    function __init__()
        @initcxx
    end
end # module Lib

using .Lib

cxxtime() = Lib.time_since_epoch()
juliatime(time = now()) = 1000 * (Dates.value(time) - Dates.UNIXEPOCH)

end


text:
# ## What is specialization?

# Example from Ref. 1.

f(x) = x * sizeof(x)
#----------------------------------------------------------------------------
@code_typed f(2)
#----------------------------------------------------------------------------
@code_typed f(2.0)

# ## References
# 1. [JuliaCon 2019 | If Runtime isn't Funtime: Controlling Compile-time Execution | Nathan Daly](https://www.youtube.com/watch?v=JCFej--XER0)


text:
abstract type AbstractModelStore end

# Required methods for subtypes:
# = initialize_storage!
# - write_result!
# - read_results
# - write_optimizer_stats!
# - read_optimizer_stats
#
# Each subtype must have a field for each instance of STORE_CONTAINERS.

function Base.empty!(store::AbstractModelStore)
    stype = typeof(store)
    for (name, type) in zip(fieldnames(stype), fieldtypes(stype))
        val = get_data_field(store, name)
        try
            empty!(val)
        catch
            @error "Base.empty! must be customized for type $stype or skipped"
            rethrow()
        end
    end
end

get_data_field(store::AbstractModelStore, type) = getfield(store, type)

function Base.isempty(store::AbstractModelStore)
    stype = typeof(store)
    for (name, type) in zip(fieldnames(stype), fieldtypes(stype))
        val = get_data_field(store, name)
        try
            !isempty(val) && return false
        catch
            @error "Base.isempty must be customized for type $stype or skipped"
            rethrow()
        end
    end

    return true
end

function list_fields(store::AbstractModelStore, container_type::Symbol)
    return keys(get_data_field(store, container_type))
end

function write_result!(store::AbstractModelStore, key, index, array)
    field = get_store_container_type(key)
    return write_result!(store, field, key, index, array)
end

function read_results(store::AbstractModelStore, key; index=nothing)
    field = get_store_container_type(key)
    return read_results(store, field, key, index=index)
end

function read_results(
    ::Type{DataFrames.DataFrame},
    store::AbstractModelStore,
    key;
    index=nothing,
)
    return read_results(store, key, index=index)
end

function list_keys(store::AbstractModelStore, container_type)
    container = get_data_field(store, container_type)
    return collect(keys(container))
end

function get_variable_value(
    store::AbstractModelStore,
    ::T,
    ::Type{U},
) where {T <: VariableType, U <: Union{PSY.Component, PSY.System}}
    return get_data_field(store, :variables)[VariableKey(T, U)]
end

function get_aux_variable_value(
    store::AbstractModelStore,
    ::T,
    ::Type{U},
) where {T <: AuxVariableType, U <: Union{PSY.Component, PSY.System}}
    return get_data_field(store, :aux_variables)[AuxVarKey(T, U)]
end

function get_dual_value(
    store::AbstractModelStore,
    ::T,
    ::Type{U},
) where {T <: ConstraintType, U <: Union{PSY.Component, PSY.System}}
    return get_data_field(store, :duals)[ConstraintKey(T, U)]
end

function get_parameter_value(
    store::AbstractModelStore,
    ::T,
    ::Type{U},
) where {T <: ParameterType, U <: Union{PSY.Component, PSY.System}}
    return get_data_field(store, :parameters)[ParameterKey(T, U)]
end


text:
using Documenter, Blue

makedocs(sitename="Blue.jl Documentation")


text:
using LearnConvection
problem = Sequential("KPP"; parameters=KPP.Parameters())

train = ["general_strat_32_profiles.jld2"]

k = 2
logγ = -1.7
distance = euclidean_distance
kernel   = get_kernel(k, logγ, 0.0, distance)

𝒟 = LearnConvection.Data.data(train, problem; D=16, N=4);
𝒢 = LearnConvection.GaussianProcess.model(𝒟; kernel = kernel)

using BenchmarkTools
@btime LearnConvection.GaussianProcess.model_output(𝒟.x[200], 𝒢)
@btime 𝒢.kernel(𝒟.x[200], 𝒢.x_train[200])

@btime 𝒢.kernel(𝒢.x_train[2], 𝒢.x_train[1])

@btime predict(𝒢, 𝒟; postprocessed=true)

get_me_true_check(𝒢, 𝒟)
get_me_true_check(𝒢, 𝒟)

using Plots
anim = animate_profile_and_model_output(𝒢, 𝒟)
gif(anim, "all.gif")

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection