text:
using ExprRules
using Test, Random

#=
TESTS NEED TO BE MORE THOROUGHLY DEFINED
=#

let
    grammar = @grammar begin
        R = I
        I = 1
    end

    @test nonterminals(grammar) == [:R,:I]
    @test return_type(grammar, 1) == :R
    @test return_type(grammar, 2) == :I
    @test child_types(grammar, 1) == [:I]
    @test child_types(grammar, 2) == Symbol[]
    @test isterminal(grammar, 1) == false
    @test isterminal(grammar, 2) == true
    @test nchildren(grammar, 1) == 1
    @test nchildren(grammar, 2) == 0
    @test max_arity(grammar) == 1

    rulenode = RuleNode(1, [RuleNode(2)])
    get_executable(rulenode, grammar)
    @test grammar[:R] == grammar.bytype[:R]
    @test grammar[:I] == grammar.bytype[:I]
    @test contains_returntype(rulenode, grammar, :R)
    @test contains_returntype(rulenode, grammar, :I)
    @test !contains_returntype(rulenode, grammar, :B)

    io = IOBuffer()
    show(io, rulenode)
    @test String(take!(io)) == "1{2}"

    display(rulenode, grammar)
end

let
    grammar = @grammar begin
        R = A
        A = B
        B = 1
    end

    rulenode = rand(RuleNode, grammar, :R)
    @test contains_returntype(rulenode, grammar, :R, 1)
    @test contains_returntype(rulenode, grammar, :R, 2)
    @test contains_returntype(rulenode, grammar, :R, 3)
    @test contains_returntype(rulenode, grammar, :A, 1) == false
    @test contains_returntype(rulenode, grammar, :A, 2)
    @test contains_returntype(rulenode, grammar, :A, 3)
    @test contains_returntype(rulenode, grammar, :B, 1) == false
    @test contains_returntype(rulenode, grammar, :B, 2) == false
    @test contains_returntype(rulenode, grammar, :B, 3)

    sample(rulenode, :R, grammar, 1)
    sample(rulenode, :R, grammar, 2)
    sample(rulenode, :R, grammar, 3)
    sample(rulenode, :A, grammar, 2)
    sample(rulenode, :A, grammar, 3)
    sample(rulenode, :B, grammar, 3)

    completed = false
    try
        sample(rulenode, :B, grammar, 2)
        completed = true
    catch
        completed = false
    end
    @test !completed #should throw an error for not found

    @test length(unique([sample(rulenode, 1) for i=1:10])) == 1
    @test length(unique([sample(rulenode, 2) for i=1:10])) <= 2
    @test length(unique([sample(rulenode, 3) for i=1:10])) <= 3

    @test length(unique([sample(NodeLoc, rulenode, 1) for i=1:10])) == 1
    @test length(unique([sample(NodeLoc, rulenode, 2) for i=1:10])) <= 2
    @test length(unique([sample(NodeLoc, rulenode, 3) for i=1:10])) <= 3

    @test sample(NodeLoc, rulenode, :R, grammar, 1) == NodeLoc(rulenode, 0)
    @test sample(NodeLoc, rulenode, :R, grammar, 2) == NodeLoc(rulenode, 0)
    @test sample(NodeLoc, rulenode, :R, grammar, 3) == NodeLoc(rulenode, 0)
    @test sample(NodeLoc, rulenode, :A, grammar, 2) == NodeLoc(rulenode, 1)
    @test sample(NodeLoc, rulenode, :A, grammar, 3) == NodeLoc(rulenode, 1)

    completed = false
    try
        sample(NodeLoc, rulenode, :A, grammar, 1)
        completed = true
    catch
        completed = false
    end
    @test !completed #should throw an error for not found
end

let
    grammar = @grammar begin
        R = R + R
        R = 1
    end

    root = RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(2)]),
                            RuleNode(1, [RuleNode(2), RuleNode(2)])])

    @test node_depth(root, root) == 1
    @test node_depth(root, root.children[1]) == 2
    @test node_depth(root, root.children[2]) == 2
    @test node_depth(root, root.children[1].children[1]) == 3
    @test node_depth(root, root.children[1].children[2]) == 3
    @test node_depth(root, root.children[2].children[1]) == 3
    @test node_depth(root, root.children[2].children[2]) == 3
end

let
    grammar = @grammar begin
        Real = x
        Real = Real * Real
        Real = f(Real)
        Real = _(rand(1.0:5.0))
        Real = A | B | g(Real)
        Real = 1 | 2 | 3
        Real = |(4:6)
        Real = |([7,8,9])
    end

    @test iseval(grammar, 1) == false
    @test iseval(grammar, 2) == false
    @test iseval(grammar, 3) == false
    @test iseval(grammar, 4) == true
    @test iseval(grammar) == true
end

x = 3
let
    Random.seed!(0)

    grammar = @grammar begin
        Real = rand()
        Real = x
        Real = Real + Real
        Real = 3.0
        Real = 1.0
        Real = sqrt(Real)
        Bool = Real > Real
        Bool = Real < Real
        Bool = true
        Real = Bool ? Real : Real
        Real = 2 / (1 + x)
        Real = 2 / (1 + Real)
    end

    @test max_arity(grammar) == 3
    @test grammar.childtypes == [
        Symbol[],
        Symbol[],
        [:Real, :Real],
        Symbol[],
        Symbol[],
        [:Real],
        [:Real, :Real],
        [:Real, :Real],
        Symbol[],
        [:Bool, :Real, :Real],
        Symbol[],
        [:Real],
    ]

    rulenode = RuleNode(3, [RuleNode(2), RuleNode(1)])
    @test length(rulenode) == 3
    @test depth(rulenode) == 2
    @test return_type(grammar, rulenode) == :Real
    @test child_types(grammar, rulenode) == [:Real,:Real]
    @test isterminal(grammar, rulenode) == false
    @test nchildren(grammar, rulenode) == 2
    @test return_type(grammar, rulenode.children[1]) == :Real
    @test child_types(grammar, rulenode.children[1]) == Symbol[]
    @test isterminal(grammar, rulenode.children[1]) == true
    @test nchildren(grammar, rulenode.children[1]) == 0
    @test  isequal(rulenode, RuleNode(3, [RuleNode(2), RuleNode(1)]))
    @test !isequal(rulenode, RuleNode(1))
    @test !isequal(rulenode, RuleNode(3, [RuleNode(1), RuleNode(2)]))
    @test string(get_executable(rulenode, grammar)) == string(:(x + rand()))

    @test contains_returntype(rulenode, grammar, :Real)
    @test !contains_returntype(rulenode, grammar, :Bool)

    rulenode = RuleNode(12, [RuleNode(3, [RuleNode(2), RuleNode(4)])])
    @test length(rulenode) == 4
    @test depth(rulenode) == 3
    @test string(get_executable(rulenode, grammar)) == string(:(2 / (1 + (x + 3.0))))

    rulenode = RuleNode(1, "potato")

    rulenode = rand(RuleNode, grammar, :Real, 3)
    hash(rulenode)

    Core.eval(rulenode, grammar)

    ex = get_executable(rulenode, grammar)
    eval(ex)
    S = SymbolTable(grammar, Main)
    interpret(S, ex)
    Core.eval(S, ex)

    Random.seed!(0)
    sample(rulenode)

    Random.seed!(1)
    sample(rulenode, :Real, grammar)

    Random.seed!(2)
    loc = sample(NodeLoc, rulenode)
    get(rulenode, loc)

    Random.seed!(3)
    insert!(rulenode, loc, rand(RuleNode, grammar, :Real, 3))

    Random.seed!(4)
    rulenode = RuleNode(3, [RuleNode(4), RuleNode(5)])
    for i in 1 : 10
        loc = sample(NodeLoc, rulenode, :Real, grammar)
    end
end

let
    grammar = @grammar begin
        R = R + R
        R = 1
        R = 2
    end

    node = RuleNode(1)
    node, worked = ExprRules._next_state!(node, grammar, 2)
    @test worked
    @test isequal(node, RuleNode(1, [RuleNode(2), RuleNode(2)]))

    node = RuleNode(1, [RuleNode(2), RuleNode(2)])
    node, worked = ExprRules._next_state!(node, grammar, 2)
    @test worked
    @test isequal(node, RuleNode(1, [RuleNode(2), RuleNode(3)]))

    node = RuleNode(1, [RuleNode(2), RuleNode(3)])
    node, worked = ExprRules._next_state!(node, grammar, 2)
    @test worked
    @test isequal(node, RuleNode(1, [RuleNode(3), RuleNode(2)]))

    node = RuleNode(1, [RuleNode(3), RuleNode(2)])
    node, worked = ExprRules._next_state!(node, grammar, 2)
    @test worked
    @test isequal(node, RuleNode(1, [RuleNode(3), RuleNode(3)]))

    node = RuleNode(1, [RuleNode(3), RuleNode(3)])
    node, worked = ExprRules._next_state!(node, grammar, 2)
    @test !worked

    ###

    node = RuleNode(1)
    for testnode in [
            RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(2)]), RuleNode(1, [RuleNode(2), RuleNode(2)])]),
            RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(2)]), RuleNode(1, [RuleNode(2), RuleNode(3)])]),
            RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(2)]), RuleNode(1, [RuleNode(3), RuleNode(2)])]),
            RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(2)]), RuleNode(1, [RuleNode(3), RuleNode(3)])]),
            RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(2)]), RuleNode(2)]),
            RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(2)]), RuleNode(3)]),
            RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(3)]), RuleNode(1, [RuleNode(2), RuleNode(2)])]),
            RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(3)]), RuleNode(1, [RuleNode(2), RuleNode(3)])]),
        ]
        node, worked = ExprRules._next_state!(node, grammar, 3)
        @test worked
        @test isequal(node, testnode)
    end

    ###
    iter = ExpressionIterator(grammar, 1, :R)
    @test all(isequal(a,b) for (a,b) in zip(collect(iter), [
        RuleNode(2), RuleNode(3)]))
    @test count_expressions(grammar, 1, :R) == 2

    iter = ExpressionIterator(grammar, 2, :R)
    x0 = iterate(iter)
    @test x0 != nothing
    item, state = x0
    @test isequal(item, RuleNode(1, [RuleNode(2), RuleNode(2)]))
    @test all(isequal(a,b) for (a,b) in zip(collect(iter), [
        RuleNode(1, [RuleNode(2), RuleNode(2)]),
        RuleNode(1, [RuleNode(2), RuleNode(3)]),
        RuleNode(1, [RuleNode(3), RuleNode(2)]),
        RuleNode(1, [RuleNode(3), RuleNode(3)]),
        RuleNode(2),
        RuleNode(3),
    ]))
    @test count_expressions(grammar, 2, :R) == 6
end

let
    grammar = @grammar begin
        R = I | F
        I = 1 | 2
        F = F + F
        F = 1.5
    end

    iter = ExpressionIterator(grammar, 2, :R)
    @test all(a == b for (a,b) in zip(collect(iter), [
        RuleNode(1, [RuleNode(3)]),
        RuleNode(1, [RuleNode(4)]),
        RuleNode(2, [RuleNode(6)]),
    ]))

    iter = ExpressionIterator(grammar, 3, :R)
    @test all(a == b for (a,b) in zip(collect(iter), [
        RuleNode(1, [RuleNode(3)]),
        RuleNode(1, [RuleNode(4)]),
        RuleNode(2, [RuleNode(5, [RuleNode(6), RuleNode(6)])]),
        RuleNode(2, [RuleNode(6)]),
    ]))

    @test count_expressions(grammar, 2, :R) == 3
    @test count_expressions(grammar, 3, :R) == 4
    @test first(iter) == RuleNode(1, [RuleNode(3)])
end

let
    grammar = @grammar begin
        R = 1
        R = 2
        R = R + R
    end

    iter = ExpressionIterator(grammar, 2, :R)
    @test all(isequal(a,b) for (a,b) in zip(collect(iter), [
        RuleNode(1),
        RuleNode(2),
        RuleNode(3, [RuleNode(1), RuleNode(1)]),
        RuleNode(3, [RuleNode(1), RuleNode(2)]),
        RuleNode(3, [RuleNode(2), RuleNode(1)]),
        RuleNode(3, [RuleNode(2), RuleNode(2)]),
    ]))

    @test count_expressions(grammar, 2, :R) == 6
    @test count_expressions(iter) == 6
    @test first(iter) == RuleNode(1)
end

let
    grammar = @grammar begin
        x = x + y
        x = y + 1
        y = x + y
        y = 1
    end

    dmap = mindepth_map(grammar)
    @test all(dmap .== [2, 1, 2, 0])
    @test mindepth(grammar, :x, dmap) == 1
    @test mindepth(grammar, :y, dmap) == 0

    node = rand(RuleNode, grammar, :y, dmap, 0)
    @test node == RuleNode(4)
    node = rand(RuleNode, grammar, :x, dmap, 1)
    @test node == RuleNode(2, [RuleNode(4)])
end

let
    grammar = @grammar begin
        b = f(x) < n
        x = true
        n = 15
    end
    node = RuleNode(1, [RuleNode(2),RuleNode(3)])
    ex = get_executable(node, grammar)
    @test ex == :(f(true) < 15)
end

let
    grammar = @grammar begin
        R = R + R
        R = 1
    end

    root = RuleNode(1, [RuleNode(1, [RuleNode(2), RuleNode(2)]),
                            RuleNode(1, [RuleNode(2), RuleNode(2)])])

    bin = NodeRecycler(20)
    recycle!(bin, root)

    RuleNode(bin, 1)
    RuleNode(bin, 1, true)
    RuleNode(bin, 1, [RuleNode(1)])
    RuleNode(bin, 1, true, [RuleNode(1)])
end

let
    grammar1 = @grammar begin
        x = x | y | a
        y = a | b
    end
    grammar2 = @grammar begin
        y = y | z
        z = a | b | c
    end

    grammar = deepcopy(grammar1)
    append!(grammar, grammar2)
end

let
    f(x1,x2,x3,x4,x5,x6,x7,x8)=x1+x2+x3+x4+x5+x6+x7+x8
    ex = :(f(1,2,3,4,5,6,7,8))
    S = SymbolTable(:f => f)
    interpret(S, ex)
    Core.eval(S, ex)
end


text:
export CTheory

struct CTheory <: AbstractTheory end

struct CTerm <: AbstractTerm
    root::Σ
    α::Union{Variable,AbstractTerm}
    β::Union{Variable,AbstractTerm}
    CTerm(root, s, t) = s >ₜ t ? new(root, t, s) : new(root, s, t)
end

function term(::CTheory, root, args)
    length(args) == 2 || throw(ArgumentError("invalid commutative term: expected 2 arguments, got $(length(args))"))
    CTerm(root, args[1], args[2])
end
Base.convert(::Type{Expr}, t::CTerm) = Expr(:call, t.root, convert(Expr, t.α), convert(Expr, t.β))

theory(::Type{CTerm}) = CTheory()
priority(::Type{CTerm}) = 50

vars(t::CTerm) = vars(t.α) ∪ vars(t.β)

Base.:(==)(a::CTerm, b::CTerm) = a.root == b.root && a.α == b.α && a.β == b.β

function (a::CTerm >ₜ b::CTerm)
    a.root == b.root || return a.root >ₑ b.root
    return a.α >ₜ b.α || (a.α == b.α && a.β >ₜ b.β)
end

Base.hash(t::CTerm, h::UInt) = hash(t.β, hash(t.α, hash(t.root, hash(CTerm, h))))

Base.map(f, t::CTerm) = CTerm(t.root, f(t.α), f(t.β))

struct CMatcher <: AbstractMatcher
    root::Σ
    s::Union{Variable,AbstractMatcher}
    t::Union{Variable,AbstractMatcher}
end

function matcher(t::CTerm, V)
    if isa(t.α, Variable) | isa(t.β, Variable) || theory(t.α) === theory(t.β)
        (cα, _) = matcher(t.α, V)
        (cβ, _) = matcher(t.β, V)
        return CMatcher(t.root, cα, cβ), V
    end

    (αβ, V1) = many_matchers([t.α, t.β], V)
    (βα, V2) = many_matchers([t.β, t.α], V)

    if length(V2) > length(V1)
        return CMatcher(t.root, βα[1], βα[2]), V2
    else
        return CMatcher(t.root, αβ[1], αβ[2]), V1
    end
end

struct CSubproblem <: AbstractSubproblem
    subproblems::Vector{Tuple{Substitution,Tuple{AbstractSubproblem,AbstractSubproblem}}}
end

function _match_c!(subproblems, σ, s, t, α, β)
    σ′ = copy(σ)

    x1 = match!(σ′, s, α)
    x1 === nothing && return

    x2 = match!(σ′, t, β)
    x2 === nothing && return

    push!(subproblems, (σ′, (x1, x2)))

    nothing
end
function match!(σ, A::CMatcher, t::CTerm)
    isa(t, CTerm) || return nothing
    A.root == t.root || return nothing

    subproblems = Tuple{Substitution,Tuple{AbstractSubproblem,AbstractSubproblem}}[]

    _match_c!(subproblems, σ, A.s, A.t, t.α, t.β)
    _match_c!(subproblems, σ, A.s, A.t, t.β, t.α)

    isempty(subproblems) && return nothing

    return CSubproblem(subproblems)
end

function compile(A::CMatcher, V)
    fn_name = gensym(:match!_c)

    subproblems = gensym(:subproblems)
    len = gensym(:len)

    cs, csf = compile(A.s, V)
    ct, ctf = compile(A.t, V)

    return fn_name, quote
        $(csf.args...)
        $(ctf.args...)
        function $fn_name(σ, t)
            isa(t, $CTerm) || return
            t.root == $(Meta.quot(A.root)) || return

            $subproblems = $(Tuple{Substitution,Tuple{AbstractSubproblem,AbstractSubproblem}})[]

            $(_compile_expr_c(subproblems, V, cs, ct))
            $(_compile_expr_c(subproblems, V, ct, cs))

            $len = $length($subproblems)
            $len == 0 && return
            $CSubproblem($subproblems)
        end
    end
end
function _compile_expr_c(subproblems, V, cs, ct)
    σ′ = gensym(:σ′)

    quote
        $σ′ = copy(σ)

        x1 = $ct($σ′, t.α)
        if x1 !== nothing
            x2 = $cs($σ′, t.β)
            if x2 !== nothing
                push!($subproblems, ($σ′, (x1, x2)))
            end
        end
    end
end

function _iterate_c_aux(p, P₀, problems, i)
    compatible(P₀, p) || return nothing
    P₁ = merge(P₀, p)
    next = _aiterate(P₁, problems...)
    next === nothing && return nothing
    (P, states) = next
    return (P, (i, P₁, states))
end
function Base.iterate(iter::Matches{CSubproblem})
    i = 1

    while i ≤ length(iter.s.subproblems)
        (P₀, problems) = iter.s.subproblems[i]
        res = _iterate_c_aux(iter.p, P₀, problems, i)
        res === nothing && (i += 1; continue)
        return res
    end

    return nothing
end

function Base.iterate(iter::Matches{CSubproblem}, (i, P₁, states))
    next = _aiterate1(P₁, iter.s.subproblems[i][2], states)

    while next === nothing
        i += 1
        i ≤ length(iter.s.subproblems) || return nothing
        (P₀, problems) = iter.s.subproblems[i]
        compatible(P₀, iter.p) || continue
        P₁ = merge(P₀, iter.p)
        next = _aiterate(P₁, problems...)
    end

    (P, states) = next
    return (P, (i, P₁, states))
end

struct CRewriter <: AbstractRewriter
    rules::Dict{Σ,Vector{Pair{CMatcher,Any}}}
end

rewriter(::CTheory) = CRewriter(Dict{Σ,Vector{Pair{CMatcher,Any}}}())
function Base.push!(rw::CRewriter, (p, b)::Pair{CTerm})
    haskey(rw.rules, p.root) || (rw.rules[p.root] = Pair{CMatcher,Any}[])
    push!(rw.rules[p.root], matcher(p) => b)
    rw
end

function compile(rw::CRewriter)
    fn_name = gensym(:rewrite_c)

    matcher_exprs = Expr[]

    rules = Dict()
    for (root, matchers) ∈ rw.rules
        rules[root] = []
        for (p, b) ∈ matchers
            fn, expr = compile(p)
            push!(matcher_exprs, expr)
            push!(rules[root], fn => b)
        end
    end

    index = gensym(:rewrite_c_map)

    fn_name, quote
        $(matcher_exprs...)
        $index = Dict($((
            :($(Meta.quot(root)) => $(:([$((:($fn => $b) for (fn, b) ∈ rs)...)])))
            for (root, rs) ∈ rules
        )...))
        function $fn_name(t)
            $haskey($rules, t.root) || return

            σ = $Substitution()

            for (match_fn, builder) ∈ $index[t.root]
                m = match_fn(σ, t)
                m === nothing && continue
                next = iterate($Matches(σ, m))
                next === nothing && continue
                σ = next[1]
                return builder(σ)::$AbstractTerm
            end

            return
        end
    end
end

ext: jl


text:
function rename(res_dir::String, new_name::String; remove_destination = false)
    mv(res_dir, joinpath(dirname(res_dir), basename(new_name)),
        remove_destination = remove_destination);
end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
# Source: https://julialang.org/blog/2016/02/iteration/#writing_multidimensional_algorithms_with_cartesianindex_iterators
# compute the "moving average" over a 3-by-3-by-... block around each element
# Origin code:
function boxcar3(A::AbstractArray)
    out = similar(A)
    R = CartesianIndices(A)
    Ifirst, Ilast = first(R), last(R)
    I1 = oneunit(Ifirst)
    for I in R
        n, s = 0, zero(eltype(out))
        for J in max(Ifirst, I-I1):min(Ilast, I+I1)
            s += A[J]
            n += 1
        end
        out[I] = s/n
    end
    out
end
#--------------------------------------------------------------
"""
    move_average(A, N)

compute the moving average over a (2N+1)-by-(2N+1)-by-... block around each element
"""
function move_average(A::AbstractArray, N::Integer)
    out = similar(A)
    R = CartesianIndices(A)
    Ifirst, Ilast = first(R), last(R)
    I1 = oneunit(Ifirst)
    for I in R
        n, s = 0, zero(eltype(out))
        for J in max(Ifirst, I-N*I1):min(Ilast, I+N*I1)
            s += A[J]
            n += 1
        end
        out[I] = s/n
    end
    out
end

max_stars_repo_path: src/moving average.jl

max_issues_repo_path: src/moving average.jl

max_forks_repo_path: src/moving average.jl


text:
## older variable visualization code below  ===========================================================

"""
    $(SIGNATURES)
Initialize empty visualizer.

DEPRECATED in favor of startDefaultVisualization
"""
function initVisualizer(;show::Bool=true)
    @warn "initVisualizer is deprecated, use startDefaultVisualization(..) instead"
    startDefaultVisualization(show=show)
end

function drawpoint!(viz,
                    sym::Symbol;
                    tf=Translation(0.0,0,0),
                    session::AbstractString="",
                    scale=0.05,
                    color=RGBA(0., 1, 0, 0.5),
                    collection::Symbol=:landmarks  )
  #

  sphere = HyperSphere(Point(0., 0, 0), scale)
  csph = GeometryData(sphere, color)
  if session == ""
    setgeometry!(viz[collection][sym], csph)
    settransform!(viz[collection][sym], tf)
  else
    sesssym=Symbol(session)
    setgeometry!(viz[sesssym][collection][sym], csph)
    settransform!(viz[sesssym][collection][sym], tf)
  end
  nothing
end

function drawpoint!(vc,
                    vert::Graphs.ExVertex,
                    topoint::Function,
                    dotwo::Bool, dothree::Bool;
                    session::AbstractString="NA"  )
  #
  den = getVertKDE(vert)
  p = Symbol(vert.label)
  pointval = topoint(den)
  if dothree
    q = convert(TransformUtils.Quaternion, Euler(pointval[4:6]...))
    drawpoint!(vc, p, tf=Translation(pointval[1:3]...), session=session)
  elseif dotwo
    drawpoint!(vc, p, tf=Translation(pointval[1],pointval[2],0.0), session=session)
  end
  nothing
end

function drawpoint!(vc,
                    vert::Graphs.ExVertex;
                    session::AbstractString="NA",
                    drawtype::Symbol=:max )
  #
  topoint = gettopoint(drawtype)
  X = getVal(vert)
  dotwo, dothree = getdotwothree(Symbol(vert.label), X)
  drawpoint!(vc, vert, topoint, dotwo, dothree, session=session)
  nothing
end

## functions to visualize a pose

function drawpose!(viz,
                   sym::Symbol;
                   tf::CTs.AbstractAffineMap=Translation(0.0,0,0)∘LinearMap(CTs.AngleAxis(0.0,0,0,1.0)),
                   session::AbstractString="",
                   collection::Symbol=:poses)
  #
  if session == ""
    setgeometry!(viz[collection][sym], Triad())
    settransform!(viz[collection][sym], tf)
  else
    sesssym=Symbol(session)
    setgeometry!(viz[sesssym][collection][sym], Triad())
    settransform!(viz[sesssym][collection][sym], tf)
  end
  nothing
end

function drawpose!(vc,
      vert::Graphs.ExVertex,
      topoint::Function,
      dotwo::Bool, dothree::Bool;
      session::AbstractString="NA"  )
  #
  den = getVertKDE(vert)
  p = Symbol(vert.label)
  pointval = topoint(den)
  tf = nothing
  if dothree
    q = convert(TransformUtils.Quaternion, Euler(pointval[4:6]...))
    tf = Translation(pointval[1:3]...)∘LinearMap(Quat(q.s,q.v...))
    drawpose!(vc, p, tf=tf, session=session)
  elseif dotwo
    tf = Translation(pointval[1],pointval[2],0.0)∘LinearMap(Rotations.AngleAxis(pointval[3],0,0,1.0))
    drawpose!(vc, p, tf=tf, session=session)
  end
  return tf
end

function drawpose!(vc,
                   vert::Graphs.ExVertex;
                   session::AbstractString="NA",
                   drawtype::Symbol=:max )
  #

  topoint = gettopoint(drawtype)
  X = getVal(vert)
  dotwo, dothree = getdotwothree(Symbol(vert.label), X)
  drawpose!(vc, vert, topoint, dotwo, dothree, session=session)
end


text:
function codegen!(ctx::Ctx)
    # Generate instructions for each function
    for decl in ctx.decls
        ctx.current_func = decl
        decl_codegen!(ctx, decl)
    end

    # Combine instructions
    asm = "; Auto generated file\n\n"

    asm *= alphalib_head()

    asm *= "\nsection .data\n"
    for inst in ctx.data
        asm *= inst * "\n"
    end

    asm *= "\nsection .text\n"
    for inst in ctx.text
        asm *= inst * "\n"
    end

    return asm
end

function decl_codegen!(ctx::Ctx, decl::Decl)
    if decl.type.kind ∈ (k_proc_t, k_fn_t)
        ctx.scratch_regs = scratchregs_new()
        ctx.torestore_regs = Set()
        ctx.code = []

        stmt_codegen!(ctx, decl.body)

        fun_name = decl.id
        frame_size = decl.nlocals * 8
        # The stack is aligned by a step of 16 bytes
        frame_size % 16 != 0 && (frame_size = floor(Int, frame_size / 16 + 1) * 16)

        # Prologue generation
        prologue = []
        pushinstr!(prologue, "global $fun_name", indent=false)
        pushinstr!(prologue, "$fun_name:", indent=false)
        pushinstr!(prologue, "push rbp")
        pushinstr!(prologue, "mov rbp, rsp")
        pushinstr!(prologue, "sub rsp, $(frame_size)")

        # Epilogue generation
        epilogue = []
        pushinstr!(epilogue, ".$fun_name.epilogue:", indent=false)
        pushinstr!(epilogue, "leave")
        pushinstr!(epilogue, "ret")

        ctx.text = vcat(ctx.text, prologue, ctx.code, epilogue)
    elseif decl.type.kind == k_int_t
        # Local variable
        exp_codegen!(ctx, decl.value)

        @assert decl.sym != nothing "The variable $(decl.id) has no symbol resolved"

        sym = sym_codegen(decl.sym)
        ctx_push!(ctx, "mov $(sym), $(regstr(decl.value.reg))")

        ctx_freescratch!(ctx, decl.value.reg)
    else
        error("Not implemented decl type $(decl.type)")
    end
end

function stmt_codegen!(ctx, stmt)
    if stmt == nothing
        return
    end

    if stmt.kind == k_stmt_exp
        exp_codegen!(ctx, stmt.exp)

        if stmt.exp.reg != nothing
            ctx_freescratch!(ctx, stmt.exp.reg)
        end
    elseif stmt.kind == k_stmt_decl
        decl_codegen!(ctx, stmt.decl)
    elseif stmt.kind == k_stmt_return
        exp = stmt.exp
        exp_codegen!(ctx, exp)

        # Move the result into rax
        ctx_push!(ctx, "mov rax, $(regstr(exp.reg))")
        ctx_push!(ctx, "jmp .$(ctx.current_func.id).epilogue")

        ctx_freescratch!(ctx, exp.reg)
    elseif stmt.kind == k_stmt_printint
        exp = stmt.exp
        exp_codegen!(ctx, exp)

        # Push preserved arg regs
        nargs = length(ctx.current_func.type.args)
        for i = 1:nargs
            ctx_push!(ctx, "push $(regstr(arg_regs[i]))")
        end

        ctx_push!(ctx, "mov rdi, $(regstr(exp.reg))")
        ctx_push!(ctx, "call alphaprintint")

        for i = nargs:-1:1
            ctx_push!(ctx, "pop $(regstr(arg_regs[i]))")
        end

        ctx_freescratch!(ctx, exp.reg)
    elseif stmt.kind == k_stmt_printstr
        # Generate data
        str = str_escape(stmt.exp)
        data = ctx_newdata!(ctx, "db $str")

        # Push preserved arg regs
        nargs = length(ctx.current_func.type.args)
        for i = 1:nargs
            ctx_push!(ctx, "push $(regstr(arg_regs[i]))")
        end

        ctx_push!(ctx, "mov rdi, $data")
        ctx_push!(ctx, "call alphaprintstr")

        for i = nargs:-1:1
            ctx_push!(ctx, "pop $(regstr(arg_regs[i]))")
        end
    elseif stmt.kind == k_stmt_scan
        # Generate data
        lsym = sym_codegen(stmt.exp.sym)

        # Push preserved arg regs
        nargs = length(ctx.current_func.type.args)
        for i = 1:nargs
            ctx_push!(ctx, "push $(regstr(arg_regs[i]))")
        end

        ctx_push!(ctx, "call alphascan")
        ctx_push!(ctx, "mov $lsym, rax")

        for i = nargs:-1:1
            ctx_push!(ctx, "pop $(regstr(arg_regs[i]))")
        end
    elseif stmt.kind == k_stmt_block
        for st in stmt.stmts
            stmt_codegen!(ctx, st)
        end
    elseif stmt.kind == k_stmt_ifelse
        labelfalse = ctx_newlabel!(ctx)
        labelend = ctx_newlabel!(ctx)

        # stmt.exp.reg is 0 if false
        exp_codegen!(ctx, stmt.exp)
        ctx_push!(ctx, "cmp $(regstr(stmt.exp.reg)), 0")
        ctx_freescratch!(ctx, stmt.exp.reg)
        ctx_push!(ctx, "je $(labelstr(labelfalse))")

        # True
        stmt_codegen!(ctx, stmt.ifbody)
        ctx_push!(ctx, "jmp $(labelstr(labelend))")

        # False
        label_codegen!(ctx, labelfalse)
        stmt_codegen!(ctx, stmt.elsebody)

        # Done
        label_codegen!(ctx, labelend)
    elseif stmt.kind == k_stmt_loop
        label_start = ctx_newlabel!(ctx)
        label_end = ctx_newlabel!(ctx)

        init, condition, iter, body = stmt.initbody, stmt.exp, stmt.iterbody, stmt.loopbody

        # Init
        stmt_codegen!(ctx, init)
        label_codegen!(ctx, label_start)

        # Test condition
        exp_codegen!(ctx, condition)
        ctx_push!(ctx, "cmp $(regstr(condition.reg)), 0")
        ctx_freescratch!(ctx, condition.reg)
        ctx_push!(ctx, "je $(labelstr(label_end))")

        # Body
        stmt_codegen!(ctx, body)

        # Iter
        stmt_codegen!(ctx, iter)
        ctx_push!(ctx, "jmp $(labelstr(label_start))")

        label_codegen!(ctx, label_end)
    else
        error("Invalid statement kind $(stmt.kind)")
    end
end

function exp_codegen!(ctx, exp)
    if exp == nothing
        return
    end

    if exp.kind == k_exp_id
        exp.reg = ctx_newscratch!(ctx)

        sym = exp.sym

        @assert sym != nothing "No symbol resolved after semantic analysis"

        ctx_push!(ctx, "mov $(regstr(exp.reg)), $(sym_codegen(sym))")
    elseif exp.kind == k_exp_int
        exp.reg = ctx_newscratch!(ctx)

        ctx_push!(ctx, "mov $(regstr(exp.reg)), $(exp.value)")
    elseif exp.kind == k_exp_neg
        l = exp.left

        exp_codegen!(ctx, l)

        # l = -l
        ctx_push!(ctx, "neg $(regstr(l.reg))")

        exp.reg = l.reg
    elseif exp.kind == k_exp_add
        l, r = exp.left, exp.right

        exp_codegen!(ctx, l)
        exp_codegen!(ctx, r)

        # l += r
        ctx_push!(ctx, "add $(regstr(l.reg)), $(regstr(r.reg))")

        ctx_freescratch!(ctx, r.reg)
        exp.reg = l.reg
    elseif exp.kind == k_exp_bool
        l, r = exp.left, exp.right

        exp_codegen!(ctx, l)
        exp_codegen!(ctx, r)

        @assert exp.operator ∈ ("and", "or", "xor", "nand", "nor") "Invalid booleand operator $(exp.operator)"

        # l op r
        ctx_push!(ctx, "$(exp.operator) $(regstr(l.reg)), $(regstr(r.reg))")

        ctx_freescratch!(ctx, r.reg)
        exp.reg = l.reg
    elseif exp.kind == k_exp_mul
        l, r = exp.left, exp.right

        exp_codegen!(ctx, l)
        exp_codegen!(ctx, r)

        # l *= r
        ctx_push!(ctx, "mov rax, $(regstr(l.reg))")
        ctx_push!(ctx, "imul $(regstr(r.reg))")
        ctx_push!(ctx, "mov $(regstr(l.reg)), rax")

        ctx_freescratch!(ctx, r.reg)
        exp.reg = l.reg
    elseif exp.kind == k_exp_div
        l, r = exp.left, exp.right

        exp_codegen!(ctx, l)
        exp_codegen!(ctx, r)

        ctx_push!(ctx, "push rdx")

        # l /= r
        ctx_push!(ctx, "mov rax, $(regstr(l.reg))")
        ctx_push!(ctx, "mov rdx, 0")
        ctx_push!(ctx, "idiv $(regstr(r.reg))")
        ctx_push!(ctx, "mov $(regstr(l.reg)), rax")

        ctx_push!(ctx, "pop rdx")

        ctx_freescratch!(ctx, r.reg)
        exp.reg = l.reg
    elseif exp.kind == k_exp_mod
        l, r = exp.left, exp.right

        exp_codegen!(ctx, l)
        exp_codegen!(ctx, r)

        ctx_push!(ctx, "push rdx")

        # l /= r
        ctx_push!(ctx, "mov rax, $(regstr(l.reg))")
        ctx_push!(ctx, "mov rdx, 0")
        ctx_push!(ctx, "idiv $(regstr(r.reg))")
        ctx_push!(ctx, "mov $(regstr(l.reg)), rdx")

        ctx_push!(ctx, "pop rdx")

        ctx_freescratch!(ctx, r.reg)
        exp.reg = l.reg
    elseif exp.kind == k_exp_set
        l, r = exp.left, exp.right

        exp_codegen!(ctx, l)
        exp_codegen!(ctx, r)

        @assert l.sym != nothing "Can't assign an rvalue"

        lsym = sym_codegen(l.sym)

        # l = r
        ctx_push!(ctx, "mov $(regstr(l.reg)), $(regstr(r.reg))")
        ctx_push!(ctx, "mov $lsym, $(regstr(l.reg))")

        ctx_freescratch!(ctx, r.reg)
        exp.reg = l.reg
    elseif exp.kind == k_exp_test
        l, r = exp.left, exp.right

        exp_codegen!(ctx, l)
        exp_codegen!(ctx, r)

        jmp_condition = test_operators[exp.operator]
        label_true = ctx_newlabel!(ctx)
        label_done = ctx_newlabel!(ctx)

        # If true, jmp to label_true
        ctx_push!(ctx, "cmp $(regstr(l.reg)), $(regstr(r.reg))")
        ctx_push!(ctx, "$(jmp_condition) $(labelstr(label_true))")

        # False
        ctx_push!(ctx, "mov $(l.reg), 0")
        ctx_push!(ctx, "jmp $(labelstr(label_done))")

        # True
        label_codegen!(ctx, label_true)
        ctx_push!(ctx, "mov $(l.reg), 1")

        # End
        label_codegen!(ctx, label_done)

        ctx_freescratch!(ctx, r.reg)
        exp.reg = l.reg
    elseif exp.kind == k_exp_call
        # Push all scratch regs
        push_scratch = scratchregs_new()
        push_scratch = [r for r in push_scratch if r ∉ ctx.scratch_regs]
        for r in push_scratch
            ctx_push!(ctx, "push $(regstr(r))")
        end

        # Push preserved arg regs
        nargs = max(length(exp.args), length(ctx.current_func.type.args))
        for i = 1:nargs
            ctx_push!(ctx, "push $(regstr(arg_regs[i]))")
        end

        for (i, arg) in enumerate(exp.args)
            exp_codegen!(ctx, arg)

            # Move the result in the arg
            argreg = arg_regs[i]
            ctx_push!(ctx, "mov $(regstr(argreg)), $(regstr(arg.reg))")

            ctx_freescratch!(ctx, arg.reg)
        end

        # Call
        ctx_push!(ctx, "call $(exp.id)")

        for i = nargs:-1:1
            ctx_push!(ctx, "pop $(regstr(arg_regs[i]))")
        end

        # Pop scratch regs
        for r in reverse(push_scratch)
            ctx_push!(ctx, "pop $(regstr(r))")
        end

        # Save return
        exp.reg = ctx_newscratch!(ctx)
        ctx_push!(ctx, "mov $(exp.reg), rax")
    else
        error("Invalid expession kind $(exp.kind)")
    end
end

function label_codegen!(ctx, label::String)
    ctx_push!(ctx, "$(label):", indent = false)
end

function label_codegen!(ctx, label::Int)
    ctx_push!(ctx, "$(labelstr(label)):", indent = false)
end

function sym_codegen(sym::Sym)
    if sym.kind == k_sym_global
        return sym.name
    elseif sym.kind == k_sym_arg
        @assert sym.position <= 5 "5 or more args are not yet supported"

        return regstr(arg_regs[sym.position])
    else
        # Local symbol
        negoffset = sym.position * 8

        return "[rbp - $negoffset]"
    end
end

# Returns the NASM compatible string
function str_escape(s::String)
    return "`$s`, 0"
end


text:
using DimensionalData, Test, Unitful, Adapt

using DimensionalData.LookupArrays, DimensionalData.Dimensions

struct CustomArray{T,N} <: AbstractArray{T,N}
    arr::Array
end

CustomArray(x::Array{T,N}) where {T,N} = CustomArray{T,N}(x)
Adapt.adapt_storage(::Type{<:CustomArray}, xs::Array) = CustomArray(xs)

Base.size(x::CustomArray, y...) = size(x.arr, y...)
Base.getindex(x::CustomArray, y...) = getindex(x.arr, y...)
Base.count(x::CustomArray) = count(x.arr)

@testset "Metadata" begin
    @test adapt(CustomArray, Metadata(:a=>"1", :b=>"2")) == NoMetadata()
end

@testset "LookupArray" begin
    l = Sampled([1:10...]; metadata=Metadata(:a=>"1", :b=>"2"))
    l1 = Adapt.adapt(CustomArray, l)
    @test parent(parent(l1)) isa CustomArray
    @test parent(parent(l1)).arr == [1:10...]
    @test metadata(l1) == NoMetadata()
    l = Categorical('a':1:'n'; metadata=Metadata(:a=>"1", :b=>"2"))
    l1 = Adapt.adapt(CustomArray, l)
    @test parent(parent(l1)) isa StepRange
    @test metadata(l1) == NoMetadata()
    l = Categorical(['a', 'e', 'n']; metadata=Metadata(:a=>"1", :b=>"2"))
    l1 = Adapt.adapt(CustomArray, l)
    @test parent(parent(l1)) isa CustomArray
    @test metadata(l1) == NoMetadata()
    l = NoLookup(Base.OneTo(10))
    l1 = Adapt.adapt(CustomArray, l)
    @test parent(parent(l1)) isa Base.OneTo
    l = AutoLookup()
    l1 = Adapt.adapt(CustomArray, l)
    @test parent(l1) isa AutoIndex
end

@testset "Dimension" begin
    d = X(Sampled([1:10...]; metadata=Metadata(:a=>"1", :b=>"2")))
    d1 = Adapt.adapt(CustomArray, d)
    @test parent(parent(d1)) isa CustomArray
    @test parent(parent(d1)).arr == [1:10...]
    @test metadata(d1) == NoMetadata()
end

@testset "DimArray" begin
    A = rand(4, 5)
    da = DimArray(A, (X, Y))
    da1 = adapt(CustomArray, da)
    @test parent(da1) isa CustomArray
    @test parent(da1).arr == A
    @test metadata(da1) == NoMetadata()
end

@testset "DimStack" begin
    A = rand(4, 5)
    B = rand(4, 5)
    ds = DimStack((a=A, b=B), (X, Y))
    ds1 = adapt(CustomArray, ds)
    @test parent(ds1[:a]) isa CustomArray
    @test parent(ds1[:a]).arr == A
    @test parent(ds1[:b]) isa CustomArray
    @test parent(ds1[:b]).arr == B
    @test metadata(ds) == NoMetadata()
end


text:
# TODO: Make tests that vary the parameters and validate the output.

plot(exp, 0, 30)
plot(sinpi, -2, 2)
plot(rand(15))
plot(rand(15, 4))

plot(x -> cos(x), 0:1)
plot(cos, 0, 1, border=false)
plot(x -> sin(x), -15:15)
plot([x -> cos(x), x -> cos(x + pi)], 0:5)
plot(x -> x^3 - 2x^2 + 3x, -5, 5)
plot(x -> tanh(x), -1.5:1.5)
plot(z -> 3z^2 - 2z^3, 0:1)

plot([1, 3, 5])
plot(1:10)
plot([1, 3, 5, 7], [13, 11, 9, 7])
plot([1:2:20], rand(10, 3))

plot(exp, 0.1:20)
logplot(exp, 0.1:20)
plot(x -> exp(-x/5), 0.01:20)
logplot(x -> exp(-x/5), 0.01:20)

plot(sinpi, -2, 2, invert=false, gridlines=false)
plot(sinpi, -2, 2, invert=true, gridlines=false)
plot(sinpi, -2, 2, invert=false, gridlines=true)
plot(sinpi, -2, 2, invert=true, gridlines=true)

# Test rounding and symoblic label display.
plot(sin, -6pi, 6pi, cols=140)


text:
using Test
using Logging
using ProgressMeter
using ValueHistories
using Plots
using Flux
using Distributions
using DistributionsAD
using ConditionalDists
using IPMeasures
using GenerativeModels
using GMExtensions

# for testing of parameter change in training
get_params(model) =  map(copy, collect(params(model)))
param_change(frozen_params, model) =
	map(x-> x[1] != x[2], zip(frozen_params, collect(params(model))))

include("saveload.jl")
include("layers.jl")
include("plotrecipes.jl")
include("callbacks.jl")
#include("reidentify.jl")
include("conv.jl")


text:
export TropicalCurve, vertices, segments, halfrays, archtrop_neighbourhoods_masked

"""
    TropicalCurve(tropicalpoly)
    TropicalCurve(tropicalpoly, newton_polygon::NewtonPolygon)

Computes a Tropical planar curve. A curce is represented as a vector of `vertices`, a vector
of `segments` between vertices and `halfrays` starting from a node in a specify direction.

You can visualize a `TropicalCurve` `curve` by using `Plots.jl` via `plot(curve)`.

    TropicalCurve(realpoly)

The polynomial will first be converted to a tropical one via [tropicalpolynomial](@ref)
and then the curve will be computed as usual.
"""
function TropicalCurve(polynomial::AbstractTropicalPolynomial, newt)
    coeffs = realcoefficients(polynomial)
    vertices = _vertices(newt, coeffs)
    segments = _segments(newt)
    halfrays = _halfrays(newt)

    TropicalCurve(polynomial, vertices, segments, halfrays)
end
TropicalCurve(p::AbstractTropicalPolynomial{<:Real}) = TropicalCurve(p, newtonpolygon(p))
TropicalCurve(p::MP.AbstractPolynomial{<:Real}) = TropicalCurve(p, tropicalpolynomial(p))

function Base.show(io::IO, curve::TropicalCurve)
    println(io, "TropicalCurve from $(curve.polynomial)")
    println(io, "vertices: ", curve.vertices)
    println(io, "segments: ", curve.segments)
    println(io, "halfrays: ", curve.halfrays)
end

function _vertices(newt::NewtonPolygon, coeffs)
    vertices = map(newt.subdivision) do simplex
        # even if our subdivision has more than 3 vertices,
        # the first 3 are sufficient to compute the node
        k, l, m = simplex
        # by duality the minimum should be attained at least three times, i.e.
        #   max(c_{i_l, j_l} + i_l x + j_l y) = max(c_{i_k, j_k} + i_k x + j_k y)
        #   max(c_{i_m, j_m} + i_m x + j_m y) = max(c_{i_k, j_k} + i_k x + j_k y)
        i_k, j_k = newt.lattices[k]
        c_k = coeffs[k]
        i_l, j_l = newt.lattices[l]
        c_l = coeffs[l]
        i_m, j_m = newt.lattices[m]
        c_m = coeffs[m]

        A = SMatrix{2,2, Float64}(i_l - i_k, i_m - i_k, j_l - j_k, j_m - j_k)
        b = SVector(c_k - c_l, c_k - c_m)
        node = inv(A) * b
    end
    vertices
end

function _segments(newt::NewtonPolygon)
    n = length(newt.subdivision)
    segments = Vector{SVector{2, Int}}()
    for i = 1:n-1, j = i+1:n
        if commonedge(newt.subdivision[i], newt.subdivision[j])
            push!(segments, SVector(i,j))
        end
    end
    segments
end

function _halfrays(newt::NewtonPolygon)
    rays = Vector{Tuple{Int64, SVector{2, Float64}}}()
    for (i, simplex) in enumerate(newt.subdivision)
        n = length(simplex)
        for k in 1:n
            a, b = simplex[k], simplex[k%n + 1]
            # we take the midpoint of a facet of the simplex and check whether it is on the
            # convex hull of the newton polygon
            classification, val = classifypoint(newt, 0.5 * newt.lattices[a] + 0.5 * newt.lattices[b])
            if classification == :hull
                push!(rays, (i, val))
            end
        end
    end
    rays
end

"""
    classifypoint(newt, p)

Check whether a points is on the convex hull of the Newton polygon, and if so return the
outer normal.
"""
function classifypoint(newt::NewtonPolygon, p::SVector{2, Float64})
    for facet in newt.facets
        i, j = facet
        a = newt.lattices[i]
        b = newt.lattices[j]
        if distance(a, b, p) < 1e-6
            # we have a point on the convex hull
            # the facets are sorted counterclockwise,
            # i.e. to get an outer normal we have to rotate 90 degree clockwise
            return (:hull, normalize(SVector(b[2] - a[2], a[1] - b[1])))
        end
    end
    (:interior, p)
end

function commonedge(a::Vector{Int}, b::Vector{Int})
    # a = sort(a)
    # b = sort(b)

    m = length(a)
    n = length(b)
    for i in 1:m
        x = a[i]
        y = a[i%m+1]
        for j = 1:n
            if (b[j] == x && b[j%n+1] == y) || (b[j%n+1] == x && b[j] == y)
                return true
            end
        end
    end
    return false
end

"""
    vertices(tropicalcurve)

The vertices of the tropical curve.
"""
vertices(c::TropicalCurve) = c.vertices

"""
    segments(tropicalcurve)

The segments of the tropical curve. Each segment is a `SVector{2}` with start and end
coordinates.
"""
segments(c::TropicalCurve) = map(s -> c.vertices[s], c.segments)

"""
    halfrays(tropicalcurve)

The halfrays of the tropical curve. Each halfray is a tuple `(coordinate, direction)`.
"""
halfrays(c::TropicalCurve) = map(pair -> (c.vertices[pair[1]], pair[2]), c.halfrays)

"""
    curve2grid(f, grid, tropicalcurve)

Map the tropical curve `tropicalcurve` on the `Grid` `grid` and call function `f` for each grid
point  on that curve.
"""
function curve2grid(f::Function, grid::Grid2D, curve::TropicalCurve)
    for segment in curve.segments
        a, b = curve.vertices[segment]
        line(f, grid, nearestpixel(grid, a), nearestpixel(grid, b))
    end

    for (raynode, raydirection) in curve.halfrays
        a = curve.vertices[raynode]
        b = a + maxlength(a, limits(grid)...) * raydirection
        line(f, grid, nearestpixel(grid, a), nearestpixel(grid, b))
    end
    nothing
end

@recipe function plot(curve::TropicalCurve)
    @assert !isempty(curve.vertices) "Curve has no vertices"

    xmin, xmax, ymin, ymax = amoeba_carcase_domain_heuristic(curve)
    if !haskey(plotattributes, :xlims)
        xlims --> (xmin, xmax)
    else
        xmin, xmax = float.(plotattributes[:xlims])
    end
    if !haskey(plotattributes, :ylims)
        ylims --> (ymin, ymax)
    else
        ymin, ymax = float.(plotattributes[:ylims])
    end

    # grid --> false
    legend --> nothing
    aspect_ratio --> :equal
    linecolor --> Colors.colorant"#CD6155"
    linewidth --> 1.0

    rays = map(curve.halfrays) do x
        raynode, raydirection = x
        a = curve.vertices[raynode]
        b, _, _ = intersection_halfray_limits(a, normalize(raydirection), xmin, xmax, ymin, ymax)
        (a, b)
    end

    for segment in curve.segments
        a, b = curve.vertices[segment]
        @series begin
            aspect_ratio --> :equal
            fill := :transparent
            seriestype := :shape
            label := ""
            [a[1], b[1]], [a[2], b[2]]
        end
    end

    for (a, b) in rays
        @series begin
            aspect_ratio --> :equal
            fill := :transparent
            seriestype := :shape
            label := ""
            [a[1], b[1]], [a[2], b[2]]
        end
    end

end


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3
reorder = 1
periods = [0,1,0]
dims = [0,0,0]
MPI.Dims_create!(nnodes, dims)
comm_cart = MPI.Cart_create(comm, dims, periods, reorder)

coords_check = Cint[-1, -1, -1]

for k = 1:dims[3], j = 1:dims[2], i = 1:dims[1]
    cart_coords = [i, j, k] .- 1  # zero-based location
    cart_rank = MPI.Cart_rank(comm_cart, cart_coords)
    @test cart_rank >= 0
    @test cart_rank < nnodes

    MPI.Cart_coords!(comm_cart, cart_rank, coords_check)
    @test all(cart_coords .== coords_check)
end

MPI.free(comm_cart)

MPI.Finalize()


text:
using MLDatasets: CIFAR10, SVHN2

loss(ŷ, y) = logitcrossentropy(ŷ, y)

num_params(model) = sum(length, Flux.params(model))

round4(x) = round(x, digits=4)

function get_data(args)
    x, y = CIFAR10.traindata()
    train_x = float(x)
    train_y = onehotbatch(y, 0:9)
    train_data = (train_x, train_y)

    test_x, test_y = CIFAR10.testdata()
    test_x = float(test_x)
    test_y = onehotbatch(test_y, 0:9)
    test_data = (test_x, test_y)

    # OOD dataset
    ood_test_x, ood_test_y = SVHN2.testdata()
    ood_test_x = float(ood_test_x)
    ood_test_y = onehotbatch(ood_test_y, 1:10)
    ood_test_data = (ood_test_x, ood_test_y)

    train_loader = DataLoader(train_data, batchsize=args.batchsize, shuffle=true, partial=false)
    test_loader = DataLoader(test_data, batchsize=args.batchsize)
    ood_test_loader = DataLoader(ood_test_data, batchsize=args.batchsize)

    return train_loader, test_loader, ood_test_loader
end


text:
using HDF5
using HDF5.Filters
using Test
using H5Zblosc, H5Zlz4, H5Zbzip2, H5Zzstd

@testset "filter" begin

# Create a new file
fn = tempname()

# Create test data
data = rand(1000, 1000)

# Open temp file for writing
f = h5open(fn, "w")

# Create datasets
dsdeflate = create_dataset(f, "deflate", datatype(data), dataspace(data),
                           chunk=(100, 100), deflate=3)

dsshufdef = create_dataset(f, "shufdef", datatype(data), dataspace(data),
                           chunk=(100, 100), shuffle=true, deflate=3)

dsfiltdef = create_dataset(f, "filtdef", datatype(data), dataspace(data),
                           chunk=(100, 100), filters=Filters.Deflate(3))

dsfiltshufdef = create_dataset(f, "filtshufdef", datatype(data), dataspace(data),
                               chunk=(100, 100), filters=[Filters.Shuffle(), Filters.Deflate(3)])

# Write data
write(dsdeflate, data)
write(dsshufdef, data)
write(dsfiltdef, data)
write(dsfiltshufdef, data)

# Test compression filters

compressionFilters = Dict(
    "blosc" => BloscFilter,
    "bzip2" => Bzip2Filter,
    "lz4" => Lz4Filter,
    "zstd" => ZstdFilter
)

for (name, filter) in compressionFilters

    ds = create_dataset(
        f, name, datatype(data), dataspace(data),
        chunk=(100,100), filters=filter()
    )
    write(ds, data)

    ds = create_dataset(
        f, "shuffle+"*name, datatype(data), dataspace(data),
        chunk=(100,100), filters=[Filters.Shuffle(), filter()]
    )
    write(ds, data)

end

# Close and re-open file for reading
close(f)
f = h5open(fn)

# Read datasets and test for equality
for name in keys(f)
    ds = f[name]
    @testset "$name" begin
        @debug "Filter Dataset" HDF5.name(ds)
        @test ds[] == data
        filters = HDF5.get_create_properties(ds).filters
        if startswith(name, "shuffle+")
            @test filters[1] isa Shuffle
            @test filters[2] isa compressionFilters[name[9:end]]
        elseif haskey(compressionFilters, name)
            @test filters[1] isa compressionFilters[name]
        end
    end
end

close(f)

end # @testset "filter"


text:
function plot_polyhedron!(region::Polyhedron, plt = current(); kwargs...)
    xs, ys = Float64[], Float64[]
    cross2 = (p1, p2) -> p1[2] * p2[1] - p1[1] * p2[2]
    points_clockwise = sort(collect(points(vrep(region))), lt = (p1, p2) -> cross2(p1, p2) <= 0)
    for point in points_clockwise
        push!(xs, point[1])
        push!(ys, point[2])
    end
    plot!(Shape(xs, ys); kwargs...)
end

function plot_piecewise_mccormick(N, umin, umax, vmin, vmax)
    DrakeVisualizer.any_open_windows() || (DrakeVisualizer.new_window(); sleep(1))
    vis = Visualizer()

    wmin = minimum(u * v for u in (umin, umax), v in (vmin, vmax))
    wmax = maximum(u * v for u in (umin, umax), v in (vmin, vmax))
    f = uvw -> uvw[1] * uvw[2] - uvw[3]
    mesh = contour_mesh(f, Vec(umin, vmin, wmin), Vec(umax, vmax, wmax))
    green = RGBA(0., 1, 0, 1.0)
    red = RGBA(1., 0, 0, 0.5)
    setgeometry!(vis[:exact], GeometryData(mesh, green))

    for i = 1 : N
        m = Model()
        @variables(m, begin
            u
            v
            w
        end)
        umin_piece = umin + (i - 1) * (umax - umin) / N
        umax_piece = umin + i * (umax - umin) / N
        envelope = mccormick_envelope(m, u, v, w, umin_piece..umax_piece, vmin..vmax)
        polyhedron_constraints(m, envelope, [u; v; w])
        poly = polyhedron(m, CDDLibrary(:exact))
        setgeometry!(vis[:mccormick][Symbol("piece$i")], GeometryData(poly, red))
    end
    vis
end

function plot_environment{T<:EnvironmentRegion}(environment::AbstractVector{T})
    plt = plot()
    for region in environment
        kwargs = if isfree(region)
            ((:lab, "free"), (:color, "blue"), (:opacity, 0.1))
        else
            ((:lab, "contact"), (:fillcolor, "black"))
        end
        plot_polyhedron!(polyhedron(region.position); kwargs...)
    end
    title!(plt, "Environment")
    plt
end

function plot_allowable_forces(region::EnvironmentRegion, id)
    plt = plot()
    plot_polyhedron!(polyhedron(region.force); opacity = 0.5, lab = "")
    title!("Allowable contact forces, region $id")
    plt
end

function plot_kinematic_regions{T<:VRepresentation}(contact_kinematic_regions::Dict{Symbol, T})
    plt = plot()
    for (name, kinematic_region) in contact_kinematic_regions
        plot_polyhedron!(polyhedron(kinematic_region); lab = string(name), opacity = 0.5)
    end
    title!("Kinematic reachability relative to CoM")
    plt
end


text:
"""
Definitions of the Pricing MDP

How show this be defined? I see two options:
 - CLASSICAL (Python?) WAY: each PMDP is given given by struct of the same type that is built
   using helper functions and stores the problem information inside the struct. (<-- What I have right now)
    - (+) I have it already this way
    - (+) I can easily have multiple instances side by side (I can have that with the other one as well)
    - (-) I am already deviating from this in HistoryReplayer, which seems it could have been a lot easier in the second option
 - MULTIPLE DISPATCH WAY: each PMDP instance is its own subtype of PMDP and to
   define a problem instance, you have to redefine abstract PMDP methods
    - (+) Seems like a natural way to write Julia
    - (+) Better decoupling of logic and instance data
    - (-) if replacing struct field (e.g. λ) with a method, the method cannot be defined by a function call (e.g. λ() = get_λ(arg) ) as that
    would call the function each time λ() is accessed. This can be addressed by putting stuff into variables and then assigning them to functions.
    - (-) sounds like work to rewrite the code this way

Conclusion: I will give it a go.
"""

"""
Values needed for defining PMDP instance
 - objective - optimization goal in the PMDP
 - edges - edges of the underlying resource graph
 - products - products are combinations of the resource edges
 - selling_period_ends -

"""

pp(m::PMDP) = m.pp

products(m::PMDP) = pp(m).P
selling_period_end(m::PMDP) = selling_period_end(pp(m))
budgets(m::PMDP) = pp(m).B
demand(m::PMDP) = pp(m).D
POMDPs.actions(m::PMDP) = pp(m).A
objective(m::PMDP) = pp(m).objective

n_resources(m::PMDP) = n_resources(pp(m))
n_products(m::PMDP) = n_products(pp(m))
n_actions(m::PMDP) = n_actions(pp(m))

empty_product(m::PMDP) = m.empty_product
empty_product_id(m::PMDP) = m.empty_product_id

POMDPs.discount(m::PMDP) = 0.99

index(m::PMDP, p::Product) = m.productindices[p]

"""
sale_prob(m::PMDP, s::State, a::Action)

Return the sale probability (Float64) of product requested in state `s` given action `a`
"""
function sale_prob end

"""
sample_customer_budget(m::PMDP, s::State, rng::AbstractRNG)

Return sampled value (of type Action) of customer budget for product requested in state `s`
"""
function sample_customer_budget end

"""
Returns user buy or no buy decision given agent selected action and user budget.
"""
user_buy(a::Action, budget::Action)::Bool = a<=budget

"""
Returns new vector of resource capacites after sale of product p
"""
# TODO: Could work inplace
reduce_capacities(c::SVector, p::Product) = c .- p

"""
Given state s, determine whether a sale of product s.p is impossible
"""
function sale_impossible(m::PMDP, s::State)::Bool
    p = product(m, s)
    s.iₚ==empty_product_id(m) || any((s.c - p) .<0.) ||  s.t >= selling_period_end(p)
end

function sale_impossible(m::PMDP, s::State, a::Action)::Bool
    a==REJECT_ACTION || sale_impossible(m, s)
end

"""
Given state s, determine whether the state is terminal in the MDP.

State is terminal if it's timestep is over the timestep limit
or if the capacity of all resources is 0.
"""
function POMDPs.isterminal(m::PMDP, s::State)::Bool
    s.t >= selling_period_end(m) || all(s.c .<= 0)
end

"""
Return an array of actions available in state s.

If product can be sold in state s, return all actions available in the MDP.
If not, return only the "impossible" action which is the first elemnt of the action array.
"""
function POMDPs.actions(m::PMDP, s::State)::AbstractArray{Action}
    actions = POMDPs.actions(m)
    if sale_impossible(m, s)
        return [actions[1]]
    else
        return actions
    end
    return actions
end

productindices(P::Array{Product{n_res}} where n_res) = Dict(zip(P, 1:length(P)))

POMDPs.initialstate(m::PMDP) = Deterministic(State(pp(m).c₀, 1, rand(demand(m)[1])))

"""
Returns the next state from given
    - state
    - action
by sampling the MDP distributions.
The most important function in the interface used by the search methods.
"""
function POMDPs.gen(m::PMDP, s::State, a::Action, rng::AbstractRNG)
    b = sample_customer_budget(m, s, rng)
    if ~sale_impossible(m, s, a) && user_buy(a, b)
        if objective(m) == :revenue
            r=a
        elseif objective(m) == :utilization
            r=sum(product(m, s))
        else
            throw(ArgumentError(string("Unknown objective: ", objective(m))))
        end
        # r = a
        c = reduce_capacities(s.c, product(m, s))
    else
        r = 0.
        c = s.c
    end
    Δt = 1
    iₚ = sample_request(m, s.t+Δt, rng)
    while iₚ==m.empty_product_id && s.t + Δt < selling_period_end(m)
        Δt += 1
        iₚ = sample_request(m, s.t+Δt, rng)
    end
    return (sp = State(c, s.t+Δt, iₚ), r = r, info=(b=b,))
end


text:
using Base: @deprecate

# BEGIN TimeZones 1.0 deprecations

@deprecate DateTime(zdt::ZonedDateTime, ::Type{Local}) DateTime(zdt)
@deprecate Date(zdt::ZonedDateTime, ::Type{Local}) Date(zdt)
@deprecate Time(zdt::ZonedDateTime, ::Type{Local}) Time(zdt)

# END TimeZones 1.0 deprecations


text:
@doc raw"""
    ProductManifold{𝔽,TM<:Tuple} <: AbstractManifold{𝔽}

Product manifold $M_1 × M_2 × …  × M_n$ with product geometry.

# Constructor

    ProductManifold(M_1, M_2, ..., M_n)

generates the product manifold $M_1 × M_2 × … × M_n$.
Alternatively, the same manifold can be contructed using the `×` operator:
`M_1 × M_2 × M_3`.
"""
struct ProductManifold{𝔽,TM<:Tuple} <: AbstractManifold{𝔽}
    manifolds::TM
end

function ProductManifold(manifolds::AbstractManifold...)
    𝔽 = ManifoldsBase._unify_number_systems((number_system.(manifolds))...)
    return ProductManifold{𝔽,typeof(manifolds)}(manifolds)
end

"""
    getindex(M::ProductManifold, i)
    M[i]

access the `i`th manifold component from the [`ProductManifold`](@ref) `M`.
"""
@inline Base.getindex(M::ProductManifold, i::Integer) = M.manifolds[i]

ProductManifold() = throw(MethodError("No method matching ProductManifold()."))

const PRODUCT_BASIS_LIST = [
    VeeOrthogonalBasis,
    DefaultBasis,
    DefaultBasis{<:Any,TangentSpaceType},
    DefaultOrthogonalBasis,
    DefaultOrthogonalBasis{<:Any,TangentSpaceType},
    DefaultOrthonormalBasis,
    DefaultOrthonormalBasis{<:Any,TangentSpaceType},
    ProjectedOrthonormalBasis{:gram_schmidt,ℝ},
    ProjectedOrthonormalBasis{:svd,ℝ},
]

"""
    ProductBasisData

A typed tuple to store tuples of data of stored/precomputed bases for a [`ProductManifold`](@ref).
"""
struct ProductBasisData{T<:Tuple}
    parts::T
end

const PRODUCT_BASIS_LIST_CACHED = [CachedBasis]

"""
    ProductMetric <: AbstractMetric

A type to represent the product of metrics for a [`ProductManifold`](@ref).
"""
struct ProductMetric <: AbstractMetric end

"""
    ProductFVectorDistribution([type::VectorBundleFibers], [x], distrs...)

Generates a random vector at point `x` from vector space (a fiber of a tangent
bundle) of type `type` using the product distribution of given distributions.

Vector space type and `x` can be automatically inferred from distributions `distrs`.
"""
struct ProductFVectorDistribution{
    TSpace<:VectorBundleFibers{<:VectorSpaceType,<:ProductManifold},
    TD<:(NTuple{N,Distribution} where {N}),
    TX,
} <: FVectorDistribution{TSpace,TX}
    type::TSpace
    x::TX
    distributions::TD
end

"""
    ProductPointDistribution(M::ProductManifold, distributions)

Product distribution on manifold `M`, combined from `distributions`.
"""
struct ProductPointDistribution{
    TM<:ProductManifold,
    TD<:(NTuple{N,Distribution} where {N}),
} <: MPointDistribution{TM}
    manifold::TM
    distributions::TD
end

"""
    ProductRetraction(retractions::AbstractRetractionMethod...)

Product retraction of `retractions`. Works on [`ProductManifold`](@ref).
"""
struct ProductRetraction{TR<:Tuple} <: AbstractRetractionMethod
    retractions::TR
end

function ProductRetraction(retractions::AbstractRetractionMethod...)
    return ProductRetraction{typeof(retractions)}(retractions)
end

"""
    InverseProductRetraction(retractions::AbstractInverseRetractionMethod...)

Product inverse retraction of `inverse retractions`. Works on [`ProductManifold`](@ref).
"""
struct InverseProductRetraction{TR<:Tuple} <: AbstractInverseRetractionMethod
    inverse_retractions::TR
end

function InverseProductRetraction(inverse_retractions::AbstractInverseRetractionMethod...)
    return InverseProductRetraction{typeof(inverse_retractions)}(inverse_retractions)
end

function allocation_promotion_function(M::ProductManifold, f, args::Tuple)
    apfs = map(MM -> allocation_promotion_function(MM, f, args), M.manifolds)
    return reduce(combine_allocation_promotion_functions, apfs)
end

"""
    ProductVectorTransport(methods::AbstractVectorTransportMethod...)

Product vector transport type of `methods`. Works on [`ProductManifold`](@ref).
"""
struct ProductVectorTransport{TR<:Tuple} <: AbstractVectorTransportMethod
    methods::TR
end

function ProductVectorTransport(methods::AbstractVectorTransportMethod...)
    return ProductVectorTransport{typeof(methods)}(methods)
end

"""
    change_representer(M::ProductManifold, ::AbstractMetric, p, X)

Since the metric on a product manifold decouples, the change of a representer can be done elementwise
"""
change_representer(::ProductManifold, ::AbstractMetric, ::Any, ::Any)

function change_representer!(M::ProductManifold, Y, G::AbstractMetric, p, X)
    map(
        (m, y, P, x) -> change_representer!(m, y, G, P, x),
        M.manifolds,
        submanifold_components(M, Y),
        submanifold_components(M, p),
        submanifold_components(M, X),
    )
    return Y
end

"""
    change_metric(M::ProductManifold, ::AbstractMetric, p, X)

Since the metric on a product manifold decouples, the change of metric can be done elementwise.
"""
change_metric(::ProductManifold, ::AbstractMetric, ::Any, ::Any)

function change_metric!(M::ProductManifold, Y, G::AbstractMetric, p, X)
    map(
        (m, y, P, x) -> change_metric!(m, y, G, P, x),
        M.manifolds,
        submanifold_components(M, Y),
        submanifold_components(M, p),
        submanifold_components(M, X),
    )
    return Y
end

"""
    check_point(M::ProductManifold, p; kwargs...)

Check whether `p` is a valid point on the [`ProductManifold`](@ref) `M`.
If `p` is not a point on `M` a [`CompositeManifoldError`](@ref) consisting of all error messages of the
components, for which the tests fail is returned.

The tolerance for the last test can be set using the `kwargs...`.
"""
function check_point(M::ProductManifold, p::Union{ProductRepr,ArrayPartition}; kwargs...)
    ts = ziptuples(Tuple(1:length(M.manifolds)), M.manifolds, submanifold_components(M, p))
    e = [(t[1], check_point(t[2:end]...; kwargs...)) for t in ts]
    errors = filter((x) -> !(x[2] === nothing), e)
    cerr = [ComponentManifoldError(er...) for er in errors]
    (length(errors) > 1) && return CompositeManifoldError(cerr)
    (length(errors) == 1) && return cerr[1]
    return nothing
end
function check_point(M::ProductManifold, p; kwargs...)
    return DomainError(
        typeof(p),
        "The point $p is not a point on $M, since currently only ProductRepr and ArrayPartition are supported types for points on arbitrary product manifolds",
    )
end

"""
    check_vector(M::ProductManifold, p, X; kwargs... )

Check whether `X` is a tangent vector to `p` on the [`ProductManifold`](@ref)
`M`, i.e. all projections to base manifolds must be respective tangent vectors.
If `X` is not a tangent vector to `p` on `M` a [`CompositeManifoldError`](@ref) consisting
of all error messages of the components, for which the tests fail is returned.

The tolerance for the last test can be set using the `kwargs...`.
"""
function check_vector(
    M::ProductManifold,
    p::Union{ProductRepr,ArrayPartition},
    X::Union{ProductRepr,ArrayPartition};
    kwargs...,
)
    ts = ziptuples(
        Tuple(1:length(M.manifolds)),
        M.manifolds,
        submanifold_components(M, p),
        submanifold_components(M, X),
    )
    e = [(t[1], check_vector(t[2:end]...; kwargs...)) for t in ts]
    errors = filter(x -> !(x[2] === nothing), e)
    cerr = [ComponentManifoldError(er...) for er in errors]
    (length(errors) > 1) && return CompositeManifoldError(cerr)
    (length(errors) == 1) && return cerr[1]
    return nothing
end
function check_vector(M::ProductManifold, p, X; kwargs...)
    return DomainError(
        typeof(X),
        "The vector $X is not a tangent vector to any tangent space on $M, since currently only ProductRepr and ArrayPartition are supported types for tangent vectors on arbitrary product manifolds",
    )
end

for TP in [ProductRepr, ArrayPartition]
    eval(
        quote
            function copyto!(M::ProductManifold, q::$TP, p::$TP)
                map(
                    copyto!,
                    M.manifolds,
                    submanifold_components(q),
                    submanifold_components(p),
                )
                return q
            end
            function copyto!(M::ProductManifold, Y::$TP, p::$TP, X::$TP)
                map(
                    copyto!,
                    M.manifolds,
                    submanifold_components(Y),
                    submanifold_components(p),
                    submanifold_components(X),
                )
                return Y
            end
        end,
    )
end

@doc raw"""
    cross(M, N)
    cross(M1, M2, M3,...)

Return the [`ProductManifold`](@ref) For two [`AbstractManifold`](@ref)s `M` and `N`,
where for the case that one of them is a [`ProductManifold`](@ref) itself,
the other is either prepended (if `N` is a product) or appenden (if `M`) is.
If both are product manifold, they are combined into one product manifold,
keeping the order.

For the case that more than one is a product manifold of these is build with the
same approach as above
"""
cross(::AbstractManifold...)
LinearAlgebra.cross(M1::AbstractManifold, M2::AbstractManifold) = ProductManifold(M1, M2)
function LinearAlgebra.cross(M1::ProductManifold, M2::AbstractManifold)
    return ProductManifold(M1.manifolds..., M2)
end
function LinearAlgebra.cross(M1::AbstractManifold, M2::ProductManifold)
    return ProductManifold(M1, M2.manifolds...)
end
function LinearAlgebra.cross(M1::ProductManifold, M2::ProductManifold)
    return ProductManifold(M1.manifolds..., M2.manifolds...)
end

@doc raw"""
    distance(M::ProductManifold, p, q)

Compute the distance between two points `p` and `q` on the [`ProductManifold`](@ref) `M`, which is
the 2-norm of the elementwise distances on the internal manifolds that build `M`.
"""
function distance(M::ProductManifold, p, q)
    return sqrt(
        sum(
            map(
                distance,
                M.manifolds,
                submanifold_components(M, p),
                submanifold_components(M, q),
            ) .^ 2,
        ),
    )
end

@doc raw"""
    exp(M::ProductManifold, p, X)

compute the exponential map from `p` in the direction of `X` on the [`ProductManifold`](@ref) `M`,
which is the elementwise exponential map on the internal manifolds that build `M`.
"""
exp(::ProductManifold, ::Any...)
function Base.exp(M::ProductManifold, p::ProductRepr, X::ProductRepr)
    return ProductRepr(
        map(
            exp,
            M.manifolds,
            submanifold_components(M, p),
            submanifold_components(M, X),
        )...,
    )
end
function Base.exp(M::ProductManifold, p::ArrayPartition, X::ArrayPartition)
    return ArrayPartition(
        map(
            exp,
            M.manifolds,
            submanifold_components(M, p),
            submanifold_components(M, X),
        )...,
    )
end

function exp!(M::ProductManifold, q, p, X)
    map(
        exp!,
        M.manifolds,
        submanifold_components(M, q),
        submanifold_components(M, p),
        submanifold_components(M, X),
    )
    return q
end

@doc raw"""
    flat(M::ProductManifold, p, X::FVector{TangentSpaceType})

use the musical isomorphism to transform the tangent vector `X` from the tangent space at
`p` on the [`ProductManifold`](@ref) `M` to a cotangent vector.
This can be done elementwise for every entry of `X` (with respect to the corresponding
entry in `p`) separately.
"""
flat(::ProductManifold, ::Any...)

function get_basis(M::ProductManifold, p, B::AbstractBasis)
    parts = map(t -> get_basis(t..., B), ziptuples(M.manifolds, submanifold_components(p)))
    return CachedBasis(B, ProductBasisData(parts))
end
function get_basis(M::ProductManifold, p, B::CachedBasis)
    return invoke(get_basis, Tuple{AbstractManifold,Any,CachedBasis}, M, p, B)
end
function get_basis(M::ProductManifold, p, B::DiagonalizingOrthonormalBasis)
    vs = map(
        ziptuples(
            M.manifolds,
            submanifold_components(p),
            submanifold_components(B.frame_direction),
        ),
    ) do t
        return get_basis(t[1], t[2], DiagonalizingOrthonormalBasis(t[3]))
    end
    return CachedBasis(B, ProductBasisData(vs))
end
for BT in PRODUCT_BASIS_LIST
    eval(quote
        @invoke_maker 3 AbstractBasis get_basis(M::ProductManifold, p, B::$BT)
    end)
end

"""
    get_component(M::ProductManifold, p, i)

Get the `i`th component of a point `p` on a [`ProductManifold`](@ref) `M`.
"""
function get_component(M::ProductManifold, p, i)
    return submanifold_component(M, p, i)
end

function get_coordinates(
    M::ProductManifold,
    p,
    X,
    B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData},
) where {𝔽}
    reps = map(
        get_coordinates,
        M.manifolds,
        submanifold_components(p),
        submanifold_components(X),
        B.data.parts,
    )
    return vcat(reps...)
end

for BT in PRODUCT_BASIS_LIST_CACHED
    eval(
        quote
            @invoke_maker 4 (
                CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData} where {𝔽}
            ) get_coordinates(M::ProductManifold, p, X, B::$BT)
        end,
    )
end
eval(
    quote
        @invoke_maker 1 AbstractManifold get_coordinates(
            M::ProductManifold,
            e::Identity,
            X,
            B::VeeOrthogonalBasis,
        )
    end,
)

function get_coordinates(M::ProductManifold, p, X, B::AbstractBasis)
    reps = map(
        t -> get_coordinates(t..., B),
        ziptuples(M.manifolds, submanifold_components(M, p), submanifold_components(M, X)),
    )
    return vcat(reps...)
end
for BT in PRODUCT_BASIS_LIST
    eval(
        quote
            @invoke_maker 4 AbstractBasis get_coordinates(M::ProductManifold, p, X, B::$BT)
        end,
    )
end

function get_coordinates!(M::ProductManifold, Xⁱ, p, X, B::AbstractBasis)
    dim = manifold_dimension(M)
    @assert length(Xⁱ) == dim
    i = one(dim)
    ts = ziptuples(M.manifolds, submanifold_components(M, p), submanifold_components(M, X))
    for t in ts
        SM = first(t)
        dim = manifold_dimension(SM)
        tXⁱ = @inbounds view(Xⁱ, i:(i + dim - 1))
        get_coordinates!(SM, tXⁱ, Base.tail(t)..., B)
        i += dim
    end
    return Xⁱ
end
function get_coordinates!(
    M::ProductManifold,
    Xⁱ,
    p,
    X,
    B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData},
) where {𝔽}
    dim = manifold_dimension(M)
    @assert length(Xⁱ) == dim
    i = one(dim)
    ts = ziptuples(
        M.manifolds,
        submanifold_components(M, p),
        submanifold_components(M, X),
        B.data.parts,
    )
    for t in ts
        SM = first(t)
        dim = manifold_dimension(SM)
        tXⁱ = @inbounds view(Xⁱ, i:(i + dim - 1))
        get_coordinates!(SM, tXⁱ, Base.tail(t)...)
        i += dim
    end
    return Xⁱ
end

for BT in PRODUCT_BASIS_LIST_CACHED
    eval(
        quote
            @invoke_maker 5 (
                CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData} where {𝔽}
            ) get_coordinates!(M::ProductManifold, Xⁱ, p, X, B::$BT)
        end,
    )
end
for BT in PRODUCT_BASIS_LIST
    eval(
        quote
            @invoke_maker 5 AbstractBasis get_coordinates!(
                M::ProductManifold,
                Xⁱ,
                p,
                X,
                B::$BT,
            )
        end,
    )
end
eval(
    quote
        @invoke_maker 1 AbstractManifold get_coordinates!(
            M::ProductManifold,
            Y,
            e::Identity,
            X,
            B::VeeOrthogonalBasis,
        )
    end,
)

function _get_dim_ranges(dims::NTuple{N,Any}) where {N}
    dims_acc = accumulate(+, vcat(1, SVector(dims)))
    return ntuple(i -> (dims_acc[i]:(dims_acc[i] + dims[i] - 1)), Val(N))
end

eval(
    quote
        @invoke_maker 1 AbstractManifold get_vector(
            M::ProductManifold,
            e::Identity,
            X,
            B::VeeOrthogonalBasis,
        )
    end,
)
for TP in [ProductRepr, ArrayPartition]
    eval(
        quote
            @invoke_maker 4 (
                CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData} where {𝔽}
            ) get_vector(
                M::ProductManifold,
                p::$TP,
                X,
                B::CachedBasis{ℝ,<:AbstractBasis{ℝ},<:ProductBasisData},
            )
            function get_vector(
                M::ProductManifold,
                p::$TP,
                X,
                B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData},
            ) where {𝔽}
                N = number_of_components(M)
                dims = map(manifold_dimension, M.manifolds)
                dim_ranges = _get_dim_ranges(dims)
                parts = ntuple(N) do i
                    return get_vector(
                        M.manifolds[i],
                        submanifold_component(p, i),
                        view(X, dim_ranges[i]),
                        B.data.parts[i],
                    )
                end
                return $TP(parts)
            end
            function get_vector(
                M::ProductManifold,
                p::$TP,
                X,
                B::AbstractBasis{𝔽,TangentSpaceType},
            ) where {𝔽}
                N = number_of_components(M)
                dims = map(manifold_dimension, M.manifolds)
                dim_ranges = _get_dim_ranges(dims)
                parts = ntuple(N) do i
                    return get_vector(
                        M.manifolds[i],
                        submanifold_component(p, i),
                        view(X, dim_ranges[i]),
                        B,
                    )
                end
                return $TP(parts)
            end
            function get_vector(M::ProductManifold, p::$TP, Xⁱ, B::VeeOrthogonalBasis)
                dim = manifold_dimension(M)
                @assert length(Xⁱ) == dim
                i = one(dim)
                ts = ziptuples(M.manifolds, submanifold_components(M, p))
                mapped = map(ts) do t
                    dim = manifold_dimension(first(t))
                    tXⁱ = @inbounds view(Xⁱ, i:(i + dim - 1))
                    i += dim
                    return get_vector(t..., tXⁱ, B)
                end
                return $TP(mapped...)
            end
        end,
    )
end

function get_vector!(M::ProductManifold, X, p, Xⁱ, B::AbstractBasis)
    dims = map(manifold_dimension, M.manifolds)
    @assert length(Xⁱ) == sum(dims)
    dim_ranges = _get_dim_ranges(dims)
    tXⁱ = map(dr -> (@inbounds view(Xⁱ, dr)), dim_ranges)
    ts = ziptuples(
        M.manifolds,
        submanifold_components(M, X),
        submanifold_components(M, p),
        tXⁱ,
    )
    map(ts) do t
        return get_vector!(t..., B)
    end
    return X
end
function get_vector!(
    M::ProductManifold,
    X,
    p,
    Xⁱ,
    B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData},
) where {𝔽}
    dims = map(manifold_dimension, M.manifolds)
    @assert length(Xⁱ) == sum(dims)
    dim_ranges = _get_dim_ranges(dims)
    tXⁱ = map(dr -> (@inbounds view(Xⁱ, dr)), dim_ranges)
    ts = ziptuples(
        M.manifolds,
        submanifold_components(M, X),
        submanifold_components(M, p),
        tXⁱ,
        B.data.parts,
    )
    map(ts) do t
        return get_vector!(t...)
    end
    return X
end
eval(
    quote
        @invoke_maker 1 AbstractManifold get_vector!(
            M::ProductManifold,
            Xⁱ,
            e::Identity,
            X,
            B::VeeOrthogonalBasis,
        )
    end,
)

for BT in PRODUCT_BASIS_LIST
    eval(
        quote
            @invoke_maker 5 AbstractBasis get_vector!(M::ProductManifold, X, p, Xⁱ, B::$BT)
        end,
    )
end
function get_vector!(M::ProductManifold, Y, p, X, B::CachedBasis)
    return error(
        "get_vector! called on $M with an incorrect CachedBasis. Expected a CachedBasis with ProductBasisData, given $B",
    )
end

function get_vectors(
    M::ProductManifold,
    p::ProductRepr,
    B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData},
) where {𝔽}
    N = number_of_components(M)
    xparts = submanifold_components(p)
    BVs = map(t -> get_vectors(t...), ziptuples(M.manifolds, xparts, B.data.parts))
    zero_tvs = map(t -> zero_vector(t...), ziptuples(M.manifolds, xparts))
    vs = typeof(ProductRepr(zero_tvs...))[]
    for i in 1:N, k in 1:length(BVs[i])
        push!(vs, ProductRepr(zero_tvs[1:(i - 1)]..., BVs[i][k], zero_tvs[(i + 1):end]...))
    end
    return vs
end
function get_vectors(
    M::ProductManifold,
    p::ArrayPartition,
    B::CachedBasis{𝔽,<:AbstractBasis{𝔽},<:ProductBasisData},
) where {𝔽}
    N = number_of_components(M)
    xparts = submanifold_components(p)
    BVs = map(t -> get_vectors(t...), ziptuples(M.manifolds, xparts, B.data.parts))
    zero_tvs = map(t -> zero_vector(t...), ziptuples(M.manifolds, xparts))
    vs = typeof(ArrayPartition(zero_tvs...))[]
    for i in 1:N, k in 1:length(BVs[i])
        push!(
            vs,
            ArrayPartition(zero_tvs[1:(i - 1)]..., BVs[i][k], zero_tvs[(i + 1):end]...),
        )
    end
    return vs
end

"""
    getindex(p, M::ProductManifold, i::Union{Integer,Colon,AbstractVector})
    p[M::ProductManifold, i]

Access the element(s) at index `i` of a point `p` on a [`ProductManifold`](@ref) `M` by
linear indexing.
See also [Array Indexing](https://docs.julialang.org/en/v1/manual/arrays/#man-array-indexing-1) in Julia.
"""
Base.@propagate_inbounds function Base.getindex(
    p::ProductRepr,
    M::ProductManifold,
    i::Union{Integer,Colon,AbstractVector,Val},
)
    return get_component(M, p, i)
end
Base.@propagate_inbounds function Base.getindex(
    p::ArrayPartition,
    M::ProductManifold,
    i::Union{Integer,Colon,AbstractVector,Val},
)
    return get_component(M, p, i)
end

@doc raw"""
    injectivity_radius(M::ProductManifold)
    injectivity_radius(M::ProductManifold, x)

Compute the injectivity radius on the [`ProductManifold`](@ref), which is the
minimum of the factor manifolds.
"""
injectivity_radius(::ProductManifold, ::Any...)
function injectivity_radius(M::ProductManifold, p)
    return min(map(injectivity_radius, M.manifolds, submanifold_components(M, p))...)
end
function injectivity_radius(M::ProductManifold, p, m::AbstractRetractionMethod)
    return min(
        map(
            (lM, lp) -> injectivity_radius(lM, lp, m),
            M.manifolds,
            submanifold_components(M, p),
        )...,
    )
end
function injectivity_radius(M::ProductManifold, p, m::ProductRetraction)
    return min(
        map(
            (lM, lp, lm) -> injectivity_radius(lM, lp, lm),
            M.manifolds,
            submanifold_components(M, p),
            m.retractions,
        )...,
    )
end
eval(
    quote
        @invoke_maker 3 AbstractRetractionMethod injectivity_radius(
            M::ProductManifold,
            p,
            B::ExponentialRetraction,
        )
    end,
)
injectivity_radius(M::ProductManifold) = min(map(injectivity_radius, M.manifolds)...)
function injectivity_radius(M::ProductManifold, m::AbstractRetractionMethod)
    return min(map(manif -> injectivity_radius(manif, m), M.manifolds)...)
end
function injectivity_radius(M::ProductManifold, m::ProductRetraction)
    return min(map((lM, lm) -> injectivity_radius(lM, lm), M.manifolds, m.retractions)...)
end
eval(
    quote
        @invoke_maker 2 AbstractRetractionMethod injectivity_radius(
            M::ProductManifold,
            B::ExponentialRetraction,
        )
    end,
)

@doc raw"""
    inner(M::ProductManifold, p, X, Y)

compute the inner product of two tangent vectors `X`, `Y` from the tangent space
at `p` on the [`ProductManifold`](@ref) `M`, which is just the sum of the
internal manifolds that build `M`.
"""
function inner(M::ProductManifold, p, X, Y)
    subproducts = map(
        inner,
        M.manifolds,
        submanifold_components(M, p),
        submanifold_components(M, X),
        submanifold_components(M, Y),
    )
    return sum(subproducts)
end

@doc raw"""
    inverse_retract(M::ProductManifold, p, q, m::InverseProductRetraction)

Compute the inverse retraction from `p` with respect to `q` on the [`ProductManifold`](@ref)
`M` using an [`InverseProductRetraction`](@ref), which by default encapsulates a inverse
retraction for each manifold of the product. Then this method is performed elementwise,
so the encapsulated inverse retraction methods have to be available per factor.
"""
inverse_retract(::ProductManifold, ::Any, ::Any, ::Any, ::InverseProductRetraction)

function inverse_retract!(M::ProductManifold, X, p, q, method::InverseProductRetraction)
    map(
        inverse_retract!,
        M.manifolds,
        submanifold_components(M, X),
        submanifold_components(M, p),
        submanifold_components(M, q),
        method.inverse_retractions,
    )
    return X
end

default_metric_dispatch(::ProductManifold, ::ProductMetric) = Val(true)

function Base.isapprox(M::ProductManifold, p, q; kwargs...)
    return all(
        t -> isapprox(t...; kwargs...),
        ziptuples(M.manifolds, submanifold_components(M, p), submanifold_components(M, q)),
    )
end
function Base.isapprox(M::ProductManifold, p, X, Y; kwargs...)
    return all(
        t -> isapprox(t...; kwargs...),
        ziptuples(
            M.manifolds,
            submanifold_components(M, p),
            submanifold_components(M, X),
            submanifold_components(M, Y),
        ),
    )
end

@doc raw"""
    log(M::ProductManifold, p, q)

Compute the logarithmic map from `p` to `q` on the [`ProductManifold`](@ref) `M`,
which can be computed using the logarithmic maps of the manifolds elementwise.
"""
log(::ProductManifold, ::Any...)

function Base.log(M::ProductManifold, p::ProductRepr, q::ProductRepr)
    return ProductRepr(
        map(
            log,
            M.manifolds,
            submanifold_components(M, p),
            submanifold_components(M, q),
        )...,
    )
end
function Base.log(M::ProductManifold, p::ArrayPartition, q::ArrayPartition)
    return ArrayPartition(
        map(
            log,
            M.manifolds,
            submanifold_components(M, p),
            submanifold_components(M, q),
        )...,
    )
end

function log!(M::ProductManifold, X, p, q)
    map(
        log!,
        M.manifolds,
        submanifold_components(M, X),
        submanifold_components(M, p),
        submanifold_components(M, q),
    )
    return X
end

@doc raw"""
    manifold_dimension(M::ProductManifold)

Return the manifold dimension of the [`ProductManifold`](@ref), which is the sum of the
manifold dimensions the product is made of.
"""
manifold_dimension(M::ProductManifold) = mapreduce(manifold_dimension, +, M.manifolds)

function mid_point!(M::ProductManifold, q, p1, p2)
    map(
        mid_point!,
        M.manifolds,
        submanifold_components(M, q),
        submanifold_components(M, p1),
        submanifold_components(M, p2),
    )
    return q
end

@doc raw"""
    norm(M::ProductManifold, p, X)

Compute the norm of `X` from the tangent space of `p` on the [`ProductManifold`](@ref),
i.e. from the element wise norms the 2-norm is computed.
"""
function LinearAlgebra.norm(M::ProductManifold, p, X)
    norms_squared = (
        map(
            norm,
            M.manifolds,
            submanifold_components(M, p),
            submanifold_components(M, X),
        ) .^ 2
    )
    return sqrt(sum(norms_squared))
end

"""
    number_of_components(M::ProductManifold{<:NTuple{N,Any}}) where {N}

Calculate the number of manifolds multiplied in the given [`ProductManifold`](@ref) `M`.
"""
number_of_components(::ProductManifold{𝔽,<:NTuple{N,Any}}) where {𝔽,N} = N

function ProductFVectorDistribution(
    type::VectorBundleFibers{<:VectorSpaceType,<:ProductManifold},
    p::Union{AbstractArray,AbstractManifoldPoint,ProductRepr},
    distributions::FVectorDistribution...,
)
    return ProductFVectorDistribution{typeof(type),typeof(distributions),typeof(p)}(
        type,
        p,
        distributions,
    )
end
function ProductFVectorDistribution(
    type::VectorBundleFibers{<:VectorSpaceType,<:ProductManifold},
    distributions::FVectorDistribution...,
)
    p = ProductRepr(map(d -> support(d).point, distributions))
    return ProductFVectorDistribution(type, p, distributions...)
end
function ProductFVectorDistribution(distributions::FVectorDistribution...)
    M = ProductManifold(map(d -> support(d).space.manifold, distributions)...)
    fiber = support(distributions[1]).space.fiber
    if !all(d -> support(d).space.fiber == fiber, distributions)
        error(
            "Not all distributions have support in vector spaces of the same type, which is currently not supported",
        )
    end
    # Probably worth considering sum spaces in the future?
    x = ProductRepr(map(d -> support(d).point, distributions)...)
    return ProductFVectorDistribution(VectorBundleFibers(fiber, M), x, distributions...)
end

function ProductPointDistribution(M::ProductManifold, distributions::MPointDistribution...)
    return ProductPointDistribution{typeof(M),typeof(distributions)}(M, distributions)
end
function ProductPointDistribution(distributions::MPointDistribution...)
    M = ProductManifold(map(d -> support(d).manifold, distributions)...)
    return ProductPointDistribution(M, distributions...)
end

function project(M::ProductManifold, p::ProductRepr)
    return ProductRepr(map(project, M.manifolds, submanifold_components(M, p))...)
end
function project(M::ProductManifold, p::ArrayPartition)
    return ArrayPartition(map(project, M.manifolds, submanifold_components(M, p))...)
end

function project!(M::ProductManifold, q, p)
    map(project!, M.manifolds, submanifold_components(M, q), submanifold_components(M, p))
    return q
end

function project(M::ProductManifold, p::ProductRepr, X::ProductRepr)
    return ProductRepr(
        map(
            project,
            M.manifolds,
            submanifold_components(M, p),
            submanifold_components(M, X),
        )...,
    )
end
function project(M::ProductManifold, p::ArrayPartition, X::ArrayPartition)
    return ArrayPartition(
        map(
            project,
            M.manifolds,
            submanifold_components(M, p),
            submanifold_components(M, X),
        )...,
    )
end

function project!(M::ProductManifold, Y, p, X)
    map(
        project!,
        M.manifolds,
        submanifold_components(M, Y),
        submanifold_components(M, p),
        submanifold_components(M, X),
    )
    return Y
end

function Random.rand(rng::AbstractRNG, d::ProductPointDistribution)
    return ProductRepr(map(d -> rand(rng, d), d.distributions)...)
end
function Random.rand(rng::AbstractRNG, d::ProductFVectorDistribution)
    return ProductRepr(map(d -> rand(rng, d), d.distributions)...)
end

function Distributions._rand!(
    rng::AbstractRNG,
    d::ProductPointDistribution,
    x::AbstractArray{<:Number},
)
    return copyto!(x, rand(rng, d))
end
function Distributions._rand!(rng::AbstractRNG, d::ProductPointDistribution, p::ProductRepr)
    map(
        t -> Distributions._rand!(rng, t[1], t[2]),
        d.distributions,
        submanifold_components(d.manifold, p),
    )
    return p
end
function Distributions._rand!(
    rng::AbstractRNG,
    d::ProductFVectorDistribution,
    v::AbstractArray{<:Number},
)
    return copyto!(v, rand(rng, d))
end
function Distributions._rand!(
    rng::AbstractRNG,
    d::ProductFVectorDistribution,
    X::ProductRepr,
)
    map(
        t -> Distributions._rand!(rng, t[1], t[2]),
        d.distributions,
        submanifold_components(d.space.manifold, X),
    )
    return X
end

@doc raw"""
    retract(M::ProductManifold, p, X, m::ProductRetraction)

Compute the retraction from `p` with tangent vector `X` on the [`ProductManifold`](@ref) `M`
using an [`ProductRetraction`](@ref), which by default encapsulates retractions of the
base manifolds. Then this method is performed elementwise, so the encapsulated retractions
method has to be one that is available on the manifolds.
"""
retract(::ProductManifold, ::Any...)

function retract!(M::ProductManifold, q, p, X, method::ProductRetraction)
    map(
        retract!,
        M.manifolds,
        submanifold_components(M, q),
        submanifold_components(M, p),
        submanifold_components(M, X),
        method.retractions,
    )
    return q
end

function representation_size(M::ProductManifold)
    return (mapreduce(m -> prod(representation_size(m)), +, M.manifolds),)
end

"""
    set_component!(M::ProductManifold, q, p, i)

Set the `i`th component of a point `q` on a [`ProductManifold`](@ref) `M` to `p`, where `p` is a point on the [`AbstractManifold`](@ref) this factor of the product manifold consists of.
"""
function set_component!(M::ProductManifold, q, p, i)
    return copyto!(submanifold_component(M, q, i), p)
end

"""
    setindex!(q, p, M::ProductManifold, i::Union{Integer,Colon,AbstractVector})
    q[M::ProductManifold,i...] = p

set the element `[i...]` of a point `q` on a [`ProductManifold`](@ref) by linear indexing to `q`.
See also [Array Indexing](https://docs.julialang.org/en/v1/manual/arrays/#man-array-indexing-1) in Julia.
"""
Base.@propagate_inbounds function Base.setindex!(
    q::Union{ProductRepr,ArrayPartition},
    p,
    M::ProductManifold,
    i::Union{Integer,Colon,AbstractVector,Val},
)
    return set_component!(M, q, p, i)
end

@doc raw"""
    sharp(M::ProductManifold, p, ξ::FVector{CotangentSpaceType})

Use the musical isomorphism to transform the cotangent vector `ξ` from the tangent space at
`p` on the [`ProductManifold`](@ref) `M` to a tangent vector.
This can be done elementwise for every entry of `ξ` (and `p`) separately
"""
sharp(::ProductManifold, ::Any...)

function _show_submanifold(io::IO, M::AbstractManifold; pre="")
    sx = sprint(show, "text/plain", M, context=io, sizehint=0)
    if occursin('\n', sx)
        sx = sprint(show, M, context=io, sizehint=0)
    end
    sx = replace(sx, '\n' => "\n$(pre)")
    print(io, pre, sx)
    return nothing
end

function _show_submanifold_range(io::IO, Ms, range; pre="")
    for i in range
        M = Ms[i]
        print(io, '\n')
        _show_submanifold(io, M; pre=pre)
    end
    return nothing
end

function _show_product_manifold_no_header(io::IO, M)
    n = length(M.manifolds)
    sz = displaysize(io)
    screen_height, screen_width = sz[1] - 4, sz[2]
    half_height = div(screen_height, 2)
    inds = 1:n
    pre = " "
    if n > screen_height
        inds = [1:half_height; (n - div(screen_height - 1, 2) + 1):n]
    end
    if n ≤ screen_height
        _show_submanifold_range(io, M.manifolds, 1:n; pre=pre)
    else
        _show_submanifold_range(io, M.manifolds, 1:half_height; pre=pre)
        print(io, "\n$(pre)⋮")
        _show_submanifold_range(
            io,
            M.manifolds,
            (n - div(screen_height - 1, 2) + 1):n;
            pre=pre,
        )
    end
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", M::ProductManifold)
    n = length(M.manifolds)
    print(io, "ProductManifold with $(n) submanifold$(n == 1 ? "" : "s"):")
    return _show_product_manifold_no_header(io, M)
end

function Base.show(io::IO, M::ProductManifold)
    return print(io, "ProductManifold(", join(M.manifolds, ", "), ")")
end

function Base.show(
    io::IO,
    mime::MIME"text/plain",
    B::CachedBasis{𝔽,T,D},
) where {𝔽,T<:AbstractBasis{𝔽},D<:ProductBasisData}
    println(io, "$(T) for a product manifold")
    for (i, cb) in enumerate(B.data.parts)
        println(io, "Basis for component $i:")
        show(io, mime, cb)
        println(io)
    end
    return nothing
end

"""
    submanifold(M::ProductManifold, i::Integer)

Extract the `i`th factor of the product manifold `M`.
"""
submanifold(M::ProductManifold, i::Integer) = M.manifolds[i]

"""
    submanifold(M::ProductManifold, i::Val)
    submanifold(M::ProductManifold, i::AbstractVector)

Extract the factor of the product manifold `M` indicated by indices in `i`.
For example, for `i` equal to `Val((1, 3))` the product manifold constructed
from the first and the third factor is returned.

The version with `AbstractVector` is not type-stable, for better preformance use `Val`.
"""
function submanifold(M::ProductManifold, i::Val)
    return ProductManifold(select_from_tuple(M.manifolds, i)...)
end
submanifold(M::ProductManifold, i::AbstractVector) = submanifold(M, Val(tuple(i...)))

Distributions.support(d::ProductPointDistribution) = MPointSupport(d.manifold)
function Distributions.support(tvd::ProductFVectorDistribution)
    return FVectorSupport(
        tvd.type,
        ProductRepr(map(d -> support(d).point, tvd.distributions)...),
    )
end

function vector_bundle_transport(::VectorSpaceType, M::ProductManifold)
    return ProductVectorTransport(map(_ -> ParallelTransport(), M.manifolds))
end
for T in [ManifoldsBase.VECTOR_TRANSPORT_DISAMBIGUATION..., AbstractVectorTransportMethod]
    eval(
        quote
            function vector_transport_direction!(M::ProductManifold, Y, p, X, d, m::$T)
                return vector_transport_direction!(
                    M,
                    Y,
                    p,
                    X,
                    d,
                    ProductVectorTransport(map(_ -> m, M.manifolds)),
                )
            end
        end,
    )
end
function vector_transport_direction!(
    M::ProductManifold,
    Y,
    p,
    X,
    d,
    m::ProductVectorTransport,
)
    map(
        vector_transport_direction!,
        M.manifolds,
        submanifold_components(M, Y),
        submanifold_components(M, p),
        submanifold_components(M, X),
        submanifold_components(M, d),
        m.methods,
    )
    return Y
end

@doc raw"""
    vector_transport_to(M::ProductManifold, p, X, q, m::ProductVectorTransport)

Compute the vector transport the tangent vector `X`at `p` to `q` on the
[`ProductManifold`](@ref) `M` using an [`ProductVectorTransport`](@ref) `m`.
This method is performed elementwise, i.e. the method `m` has to be implemented on the
base manifold.
"""
vector_transport_to(::ProductManifold, ::Any, ::Any, ::Any, ::ProductVectorTransport)

for T in [ManifoldsBase.VECTOR_TRANSPORT_DISAMBIGUATION..., AbstractVectorTransportMethod]
    eval(
        quote
            function vector_transport_to!(M::ProductManifold, Y, p, X, q, m::$T)
                return vector_transport_to!(
                    M,
                    Y,
                    p,
                    X,
                    q,
                    ProductVectorTransport(map(_ -> m, M.manifolds)),
                )
            end
        end,
    )
end
function vector_transport_to!(M::ProductManifold, Y, p, X, q, m::ProductVectorTransport)
    map(
        vector_transport_to!,
        M.manifolds,
        submanifold_components(M, Y),
        submanifold_components(M, p),
        submanifold_components(M, X),
        submanifold_components(M, q),
        m.methods,
    )
    return Y
end

function zero_vector!(M::ProductManifold, X, p)
    map(
        zero_vector!,
        M.manifolds,
        submanifold_components(M, X),
        submanifold_components(M, p),
    )
    return X
end


text:
module PushshiftRedditDistiller

using DataDeps

export RedditDataFilter, distill

include("deps.jl")
include("types.jl")
include("distiller.jl")
include("filter.jl")
include("utils.jl")

function __init__()
    init_deps()
end

end


text:
using ..Mesh.Grids: _x1, _x2, _x3
"""
    InitStateBC

Set the value at the boundary to match the `init_state_prognostic!` function. This is
mainly useful for cases where the problem has an explicit solution.

# TODO: This should be fixed later once BCs are figured out (likely want
# different things here?)
"""
struct InitStateBC end
function boundary_state!(
    ::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc::InitStateBC,
    m::AtmosModel,
    state⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    # Put cood in a NamedTuple to mimmic LocalGeometry
    init_state_prognostic!(m, state⁺, aux⁺, (coord = aux⁺.coord,), t)
end

function boundary_state!(
    ::NumericalFluxSecondOrder,
    bc::InitStateBC,
    m::AtmosModel,
    state⁺::Vars,
    diff⁺::Vars,
    hyperdiff⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    diff⁻::Vars,
    hyperdiff⁻::Vars,
    aux⁻::Vars,
    t,
    args...,
)
    # Put coord in a NamedTuple to mimmic LocalGeometry
    init_state_prognostic!(m, state⁺, aux⁺, (coord = aux⁺.coord,), t)
end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

using SlimOptim, Test

@testset "SPG" begin
    include("test_spg.jl")
end

@testset "PQN" begin
    include("test_pqn.jl")
end

@testset "Linearized bregman" begin
    include("test_bregman.jl")
end

@testset "Rosenbrock accuracy" begin
    include("test_rosen.jl")
end