text:
module Spikes

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict()

const entities = Dict{String, Type}(
    "up" => Maple.SpikesUp,
    "down" => Maple.SpikesDown,
    "left" => Maple.SpikesLeft,
    "right" => Maple.SpikesRight,
)

const triggerEntities = Dict{String, Type}(
    "up" => Maple.TriggerSpikesUp,
    "down" => Maple.TriggerSpikesDown,
    "left" => Maple.TriggerSpikesLeft,
    "right" => Maple.TriggerSpikesRight,
)

const triggerEntitiesOrig = Dict{String, Type}(
    "up" => Maple.TriggerSpikesOriginalUp,
    "down" => Maple.TriggerSpikesOriginalDown,
    "left" => Maple.TriggerSpikesOriginalLeft,
    "right" => Maple.TriggerSpikesOriginalRight,
)

const normalSpikesUnion = Union{Maple.SpikesUp, Maple.SpikesDown, Maple.SpikesLeft, Maple.SpikesRight}
const triggerSpikesUnion = Union{Maple.TriggerSpikesUp, Maple.TriggerSpikesDown, Maple.TriggerSpikesLeft, Maple.TriggerSpikesRight}
const triggerSpikesOrigUnion = Union{Maple.TriggerSpikesOriginalUp, Maple.TriggerSpikesOriginalDown, Maple.TriggerSpikesOriginalLeft, Maple.TriggerSpikesOriginalRight}
const spikesUnion = Union{normalSpikesUnion, triggerSpikesUnion, triggerSpikesOrigUnion}

const triggerSpikeColors = [
    (242, 90, 16, 255) ./ 255,
    (255, 0, 0, 255) ./ 255,
    (242, 16, 103, 255) ./ 255
]

for variant in Maple.spike_types
    for (dir, entity) in entities
        key = "Spikes ($(uppercasefirst(dir)), $(uppercasefirst(variant)))"
        placements[key] = Ahorn.EntityPlacement(
            entity,
            "rectangle",
            Dict{String, Any}(
                "type" => variant
            )
        )
    end

    if variant != "tentacles"
        for (dir, entity) in triggerEntitiesOrig
            key = "Trigger Spikes ($(uppercasefirst(dir)), $(uppercasefirst(variant)))"
            placements[key] = Ahorn.EntityPlacement(
                entity,
                "rectangle",
                Dict{String, Any}(
                    "type" => variant
                )
            )
        end
    end
end

for (dir, entity) in triggerEntities
    key = "Trigger Spikes ($(uppercasefirst(dir)), Dust)"
    placements[key] = Ahorn.EntityPlacement(
        entity,
        "rectangle"
    )
end

Ahorn.editingOptions(entity::normalSpikesUnion) = Dict{String, Any}(
    "type" => Maple.spike_types
)

Ahorn.editingOptions(entity::triggerSpikesOrigUnion) = Dict{String, Any}(
    "type" => String[variant for variant in Maple.spike_types if variant != "tentacles"]
)

const directions = Dict{String, String}(
    "spikesUp" => "up",
    "spikesDown" => "down",
    "spikesLeft" => "left",
    "spikesRight" => "right",

    "triggerSpikesUp" => "up",
    "triggerSpikesDown" => "down",
    "triggerSpikesLeft" => "left",
    "triggerSpikesRight" => "right",

    "triggerSpikesOriginalUp" => "up",
    "triggerSpikesOriginalDown" => "down",
    "triggerSpikesOriginalLeft" => "left",
    "triggerSpikesOriginalRight" => "right",
)

const offsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (4, -4),
    "down" => (4, 4),
    "left" => (-4, 4),
    "right" => (4, 4),
)

const triggerOriginalOffsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (0, 5),
    "down" => (0, -4),
    "left" => (5, 0),
    "right" => (-4, 0),
)

const rotations = Dict{String, Number}(
    "up" => 0,
    "right" => pi / 2,
    "down" => pi,
    "left" => pi * 3 / 2
)

const rotationOffsets = Dict{String, Tuple{Number, Number}}(
    "up" => (0.5, 0.25),
    "right" => (1, 0.675),
    "down" => (1.5, 1.125),
    "left" => (0, 1.675)
)

const resizeDirections = Dict{String, Tuple{Bool, Bool}}(
    "up" => (true, false),
    "down" => (true, false),
    "left" => (false, true),
    "right" => (false, true),
)

const tentacleSelectionOffsets = Dict{String, Tuple{Number, Number}}(
    "up" => (0, -8),
    "down" => (0, -8),
    "left" => (-8, 0),
    "right" => (-8, 0)
)

const spikeNames = String["spikesDown", "spikesLeft", "spikesRight", "spikesUp"]

const triggerNames = String["triggerSpikesDown", "triggerSpikesLeft", "triggerSpikesRight", "triggerSpikesUp"]
const triggerRotationOffsets = Dict{String, Tuple{Number, Number}}(
    "up" => (3, -1),
    "right" => (4, 3),
    "down" => (5, 5),
    "left" => (-1, 4),
)

const triggerOriginalNames = String["triggerSpikesOriginalDown", "triggerSpikesOriginalLeft", "triggerSpikesOriginalRight", "triggerSpikesOriginalUp"]

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::spikesUnion)
    direction = get(directions, entity.name, "up")
    theta = rotations[direction] - pi / 2

    width = Int(get(entity.data, "width", 0))
    height = Int(get(entity.data, "height", 0))

    x, y = Ahorn.position(entity)
    cx, cy = x + floor(Int, width / 2) - 8 * (direction == "left"), y + floor(Int, height / 2) - 8 * (direction == "up")

    Ahorn.drawArrow(ctx, cx, cy, cx + cos(theta) * 24, cy + sin(theta) * 24, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.selection(entity::spikesUnion)
    if haskey(directions, entity.name)
        x, y = Ahorn.position(entity)

        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))

        direction = get(directions, entity.name, "up")
        variant = get(entity.data, "type", "default")

        if variant == "tentacles"
            ox, oy = tentacleSelectionOffsets[direction]

            width = Int(get(entity.data, "width", 16))
            height = Int(get(entity.data, "height", 16))

            return Ahorn.Rectangle(x + ox, y + oy, width, height)

        else
            width = Int(get(entity.data, "width", 8))
            height = Int(get(entity.data, "height", 8))

            ox, oy = offsets[direction]

            return Ahorn.Rectangle(x + ox - 4, y + oy - 4, width, height)
        end
    end
end

function Ahorn.minimumSize(entity::spikesUnion)
    if haskey(directions, entity.name)
        variant = get(entity.data, "type", "default")

        return variant == "tentacles" ? (16, 16) : (8, 8)
    end
end

function Ahorn.resizable(entity::spikesUnion)
    if haskey(directions, entity.name)
        variant = get(entity.data, "type", "default")
        direction = get(directions, entity.name, "up")

        return resizeDirections[direction]
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::spikesUnion)
    if haskey(directions, entity.name)
        variant = entity.name in triggerNames ? "trigger" : get(entity.data, "type", "default")
        direction = get(directions, entity.name, "up")
        triggerOriginalOffset = entity.name in triggerOriginalNames ? triggerOriginalOffsets[direction] : (0, 0)

        if variant == "tentacles"
            width = get(entity.data, "width", 16)
            height = get(entity.data, "height", 16)

            for ox in 0:16:width - 16, oy in 0:16:height - 16
                drawX = ox + 16 * rotationOffsets[direction][1] + triggerOriginalOffset[1]
                drawY = oy + 16 * rotationOffsets[direction][2] + triggerOriginalOffset[2]

                Ahorn.drawSprite(ctx, "danger/tentacles00", drawX, drawY, rot=rotations[direction])
            end

            if width / 8 % 2 == 1 || height / 8 % 2 == 1
                drawX = width - 16 + 16 * rotationOffsets[direction][1] + triggerOriginalOffset[1]
                drawY = height - 16 + 16 * rotationOffsets[direction][2] + triggerOriginalOffset[2]

                Ahorn.drawSprite(ctx, "danger/tentacles00", drawX, drawY, rot=rotations[direction])
            end

        elseif variant == "trigger"
            rng = Ahorn.getSimpleEntityRng(entity)

            width = get(entity.data, "width", 8)
            height = get(entity.data, "height", 8)

            updown = direction == "up" || direction == "down"

            for ox in 0:8:width - 8, oy in 0:8:height - 8
                color1 = rand(rng, triggerSpikeColors)
                color2 = rand(rng, triggerSpikeColors)

                drawX = ox + triggerRotationOffsets[direction][1] + triggerOriginalOffset[1]
                drawY = oy + triggerRotationOffsets[direction][2] + triggerOriginalOffset[2]

                Ahorn.drawSprite(ctx, "danger/triggertentacle/wiggle_v06", drawX, drawY, rot=rotations[direction], tint=color1)
                Ahorn.drawSprite(ctx, "danger/triggertentacle/wiggle_v03", drawX + 3 * updown, drawY + 3 * !updown, rot=rotations[direction], tint=color2)
            end

        else
            width = get(entity.data, "width", 8)
            height = get(entity.data, "height", 8)

            for ox in 0:8:width - 8, oy in 0:8:height - 8
                drawX = ox + offsets[direction][1] + triggerOriginalOffset[1]
                drawY = oy + offsets[direction][2] + triggerOriginalOffset[2]

                Ahorn.drawSprite(ctx, "danger/spikes/$(variant)_$(direction)00", drawX, drawY)
            end
        end
    end
end

for (to, from) in [(Maple.SpikesUp, Maple.SpikesDown), (Maple.TriggerSpikesUp, Maple.TriggerSpikesDown), (Maple.TriggerSpikesOriginalUp, Maple.TriggerSpikesOriginalDown)]
    @eval function Ahorn.flipped(entity::$to, horizontal::Bool)
        if !horizontal
            return $from(entity.x, entity.y, entity.width, entity.type)
        end
    end

    @eval function Ahorn.flipped(entity::$from, horizontal::Bool)
        if !horizontal
            return $to(entity.x, entity.y, entity.width, entity.type)
        end
    end
end

for (to, from) in [(Maple.SpikesLeft, Maple.SpikesRight), (Maple.TriggerSpikesLeft, Maple.TriggerSpikesRight), (Maple.TriggerSpikesOriginalLeft, Maple.TriggerSpikesOriginalRight)]
    @eval function Ahorn.flipped(entity::$to, horizontal::Bool)
        if horizontal
            return $from(entity.x, entity.y, entity.height, entity.type)
        end
    end

    @eval function Ahorn.flipped(entity::$from, horizontal::Bool)
        if horizontal
            return $to(entity.x, entity.y, entity.height, entity.type)
        end
    end
end

# TODO - Rotations might need offsets

const spikesUp = [Maple.SpikesUp, Maple.TriggerSpikesUp, Maple.TriggerSpikesOriginalUp]
const spikesRight = [Maple.SpikesRight, Maple.TriggerSpikesRight, Maple.TriggerSpikesOriginalRight]
const spikesDown = [Maple.SpikesDown, Maple.TriggerSpikesDown, Maple.TriggerSpikesOriginalDown]
const spikesLeft = [Maple.SpikesLeft, Maple.TriggerSpikesLeft, Maple.TriggerSpikesOriginalLeft]

for (left, normal, right) in zip(spikesLeft, spikesUp, spikesRight)
    @eval function Ahorn.rotated(entity::$normal, steps::Int)
        if steps > 0
            return Ahorn.rotated($right(entity.x, entity.y, entity.width, entity.type), steps - 1)

        elseif steps < 0
            return Ahorn.rotated($left(entity.x, entity.y, entity.width, entity.type), steps + 1)
        end

        return entity
    end
end

for (left, normal, right) in zip(spikesUp, spikesRight, spikesDown)
    @eval function Ahorn.rotated(entity::$normal, steps::Int)
        if steps > 0
            return Ahorn.rotated($right(entity.x, entity.y, entity.height, entity.type), steps - 1)

        elseif steps < 0
            return Ahorn.rotated($left(entity.x, entity.y, entity.height, entity.type), steps + 1)
        end

        return entity
    end
end

for (left, normal, right) in zip(spikesRight, spikesDown, spikesLeft)
    @eval function Ahorn.rotated(entity::$normal, steps::Int)
        if steps > 0
            return Ahorn.rotated($right(entity.x, entity.y, entity.width, entity.type), steps - 1)

        elseif steps < 0
            return Ahorn.rotated($left(entity.x, entity.y, entity.width, entity.type), steps + 1)
        end

        return entity
    end
end

for (left, normal, right) in zip(spikesDown, spikesLeft, spikesUp)
    @eval function Ahorn.rotated(entity::$normal, steps::Int)
        if steps > 0
            return Ahorn.rotated($right(entity.x, entity.y, entity.height, entity.type), steps - 1)

        elseif steps < 0
            return Ahorn.rotated($left(entity.x, entity.y, entity.height, entity.type), steps + 1)
        end

        return entity
    end
end

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

function subscript(i)
    map(repr(i)) do c
        c   ==  '1' ? '\u2081' :
        c   ==  '2' ? '\u2082' :
        c   ==  '3' ? '\u2083' :
        c   ==  '4' ? '\u2084' :
        c   ==  '5' ? '\u2085' :
        c   ==  '6' ? '\u2086' :
        c   ==  '7' ? '\u2087' :
        c   ==  '8' ? '\u2088' :
        c   ==  '9' ? '\u2089' :
        c   ==  '0' ? '\u2080' :
        error("Unexpected character")
    end
end

function pretty_print_constraint_equation(a, b, g; char1="u", char2="λ")
    s = ""
    for (i, c) in enumerate(find(a))
        if (i != 1)
            s = s * ((a[c] < 0) ? " - " : " + ")
        end
        s = s*"$(round(abs(a[c]), 3))*$char1$(subscript(c))"
    end
    for (i, c) in enumerate(find(b))
        if (length(find(a)) != 0) || (i != 1)
            s = s * ((b[c] < 0) ? " - " : " + ")
        end
        s = s*"$(round(abs(b[c]), 3))*$char2$(subscript(c))"
    end
    s = s * " = $(round(g, 3))"
    s
end

function pretty_print_C1_row(r)
    s = "⋯ "
    for a in find(r)
        c = r[a] < 0 ? "-" : "+"
        s = s * "$c $(round(abs(r[a]), 3))*λ$(subscript(a)) "
    end
    s = s * "⋯"
    return s
end

function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2, D_, D, g_, g; show_info=false)
    # old, new, old, new...

#= herzian contact with symmetry boundary
# INFO: SUMMARY for node id 555 with dofs 1109, 1110:
# INFO: ----- Current constraint -----
# INFO: lambda coefficients in C1 matrix are:
# INFO: dof 1109: ⋯ + 0.15*λ₁₁₀₉ ⋯
# INFO: rows in constraint matrix C2 & D
# INFO: dof 1109: 0.15*u₁₁₀₉ = -0.0 <-- overconstrained dof
# INFO: ----- New constraint -----
# INFO: lambda coefficients in C1 matrix are:
# INFO: dof 1109: ⋯ + 0.15*λ₁₁₀₉ ⋯
# INFO: dof 1110: ⋯ + 0.15*λ₁₁₁₀ ⋯
# INFO: rows in constraint matrix C2 & D
# INFO: dof 1109: 0.0*u₁₅₃ - 0.0*u₁₅₄ - 0.0*u₁₅₅ + 0.15*u₁₅₆ + 0.0*u₁₁₀₉ - 0.15*u₁₁₁₀ = -0.0 <-- overconstrained dof
# INFO: dof 1110: 0.15*λ₁₁₀₉ + 0.0*λ₁₁₁₀ = -0.0
# INFO: ----- Related equations -----
# INFO: dof 1110: 0.15*λ₁₁₀₉ + 0.0*λ₁₁₁₀ = -0.0
# INFO: dof 155: 0.165*u₁₅₅ = 0.0
# INFO: algorithm 1 solved issue? false
# INFO: algorithm 2 solved issue? true
# INFO: fixed: new setting is
# INFO: dof 1109: 0.0*u₁₅₃ - 0.0*u₁₅₄ - 0.0*u₁₅₅ + 0.15*u₁₅₆ + 0.0*u₁₁₀₉ - 0.15*u₁₁₁₀ = -0.0
=#
#=
    if 555 in nodes
        info("overconstraint DIRTY HACK")
        # It is possible to selectively remove mortar constraints and the associated
        # Lagrange multiplier components in certain axis directions and replace them
        # with the Dirichlet (symmetry) conditions.
        # old configuration is dirichlet symmetry condition
        # new configuration is mortar constraint
        # 1. remove mortar constrains and associated Lagrange multiplier components
        # in dof 1109, that is, x direction of node 555.

        # works quite well
#       C1_[1109,:] = 0
#       C2_[1110,:] = C2_[1109,:]
#       C2[1110,:] = 0
#       D[1110,:] = 0

        #C1_[1110,:] = C1_[1109,:]
        C2_[1110,:] = C2_[1109,:]
        #C1_[1109,:] = 0
        C2_[1109,:] = 0

        #C1[1110,:] = 0
        #C2[1110,:] = 0
        #C2[:,1109] = 0
        D[1110,:] = 0
        #D[:,1109] = 0
        g[1110,:] = 0
        #D[:,1109] = 0
        #C2[:,1109] = 0
        #C1[1109,:] = 0
        #C1_[1110,:] = 0
        return
    end
=#

    """ Return all other dofs which connects to overconstrained dofs. """
    function get_related_dofs(dofs)
        dofs_ = Set(dofs)
        for dof in copy(dofs_)
            c = find(C2[dof,:])
            length(c) != 0 && push!(dofs_, c...)
            c = find(C2_[dof,:])
            length(c) != 0 && push!(dofs_, c...)
        end
        for dof in copy(dofs_)
            c = find(C2[dof,:])
            length(c) != 0 && push!(dofs_, c...)
            c = find(C2_[dof,:])
            length(c) != 0 && push!(dofs_, c...)
        end
        for dof in copy(dofs_)
            c = find(C2[dof,:])
            length(c) != 0 && push!(dofs_, c...)
            c = find(C2_[dof,:])
            length(c) != 0 && push!(dofs_, c...)
        end
        dofs_ = sort(collect(dofs_))
        return dofs_
    end

    """ Return true if dofs has lagrange coefficients, i.e. D is nonzero. """
    function has_lagrange_coefficients(dof::Int)
        return (countnz(D[dof, :]) != 0) || (countnz(D_[dof, :]) != 0)
    end
    function has_lagrange_coefficients(dofs::Vector{Int})
        return map(has_lagrange_coefficients, dofs)
    end

    """ Test is dof single point constraint. """
    function is_spc(C, dof)
        return countnz(C[dof, :]) == 1
    end
    function is_spc(dof::Int)
        return is_spc(C2, dof) && is_spc(C2_, dof)
    end
    function is_spc(dofs::Vector{Int})
        return map(is_spc, dofs)
    end

    function has_anything(dof::Int)
        countnz(C1[dof,:]) != 0 && return true
        countnz(C2[dof,:]) != 0 && return true
        countnz(D[dof,:]) != 0 && return true
        countnz(g[dof,:]) != 0 && return true
        return false
    end

    """ Algorithm 1. Calculate rank of overdetermined system and do LSQ if
        rank(C) equals to number of unique dofs.
    """
    function action1(node_id, dofs)
        dofs_ = get_related_dofs(intersect(dofs, all_dofs))
        # this will fail with dofs > 2 for some yet unknown reason
        length(dofs_) > 2 && return dofs_, false

        any(has_lagrange_coefficients(dofs_)) && return dofs_, false
        C = full([C2[dofs_, :]; C2_[dofs_, :]])
        d = full([g[dofs_]; g_[dofs_]])
        info("rank = $(rank(C)), dofs = $(length(dofs_))")
        rank(C) != length(dofs_) && return dofs_, false
        C2_[dofs_,:] = C2[dofs_,:] = 0
        g_[dofs_,:] = g[dofs_,:] = 0
        x = C \ d
        C2[dofs_, dofs_] = eye(length(dofs_))
        g[dofs_] = x
        return dofs_, true
    end

    function action2(node_id, dofs)
        """ If no coefficients on matrix D we can set essential boundary condition
            only on master side and impose bc weakly on slave side, i.e., remove SPC
        """
        dofs_ = intersect(dofs, all_dofs)
        length(dofs_) == 1 || return dofs_, false
        any(has_lagrange_coefficients(dofs_)) && return dofs_, false
        if is_spc(C2, dofs_) && !is_spc(C2_, dofs_)
            C1[dofs_,:] = C2[dofs_,:] = g[dofs_,:] = 0
            return dofs_, true
        elseif is_spc(C2_, dofs_) && !is_spc(C2, dofs_)
            C1_[dofs_,:] = C2_[dofs_,:] = g_[dofs_,:] = 0
            return dofs_, true
        else
            return dofs_, false
        end
    end

    function action3(node_id, dofs)
        """ Another option is to obey single point constraints
        """
        dofs_ = intersect(dofs, all_dofs)
        any(has_lagrange_coefficients(dofs_)) && return dofs_, false
        if !is_spc(C2, dofs_) && is_spc(C2_, dofs_)
            C1[dofs_,:] = C2[dofs_,:] = g[dofs_,:] = 0
            return dofs_, true
        elseif !is_spc(C2_, dofs_) && is_spc(C2, dofs_)
            C1_[dofs_,:] = C2_[dofs_,:] = g_[dofs_,:] = 0
            return dofs_, true
        else
            return dofs_, false
        end
    end

    function action4(node_id, dofs)
        """ If symmetry line, one possibility is to apply both conditions and
        eliminate lagrange multiplier. """
        dofs_ = intersect(dofs, all_dofs)
        length(dofs_) != 1 && return dofs_, false
        related_dofs = get_related_dofs(dofs_)
        for j in related_dofs
            has_anything(j) && continue
            # copy one constaint to this dof
            C1[j,:] = C1[dofs_,:]
            C2[j,:] = C2[dofs_,:]
            D[j,:] = D[dofs_,:]
            g[j,:] = g[dofs_,:]
            # make room for new constraint
            C1[dofs_,:] = 0
            C2[dofs_,:] = 0
            D[dofs_,:] = 0
            g[dofs_,:] = 0
            dofs_ = [dofs_; j]
            break
        end
        for j in related_dofs
            has_anything(j) && continue
            # set lagrange multiplier to 1
            D[j,dofs_[1]] = 1.0
            C1[j,dofs_[1]] = 1.0
            dofs_ = [dofs_; j]
            break
        end
        return dofs_, true
    end

    actions = [action1, action2, action3, action4]

    function show_lambda_coefficients(dofs, C1)
        for dof in dofs
            status = dof in all_dofs ? " <-- overconstrained dof" : ""
            r = C1[:, dof]
            r[abs(r) .< 1.0e-9] = 0
            length(nonzeros(r)) != 0 || continue
            info("dof $dof: "*pretty_print_C1_row(r))
        end
    end

    function show_rows_in_constraint_matrix(dofs, C2, D; show_status=true)
        for dof in dofs
            status = (dof in all_dofs) && show_status ? " <-- overconstrained dof" : ""
            a = C2[dof,:]
            b = D[dof,:]
            a[abs(a) .< 1.0e-12] = 0
            b[abs(b) .< 1.0e-12] = 0
            if (length(nonzeros(a)) == 0) && (length(nonzeros(b)) == 0)
                continue
            end
            info("dof $dof: "*pretty_print_constraint_equation(a, b, g[dof])*status)
        end
    end

    function show_related_equations(dofs, C2, C2_, D, D_)
        related_dofs = Set{Int64}()
        G = [C2 D]
        G_ = [C2_ D_]
        for dof in dofs
            dof in all_dofs || continue
            c = find(G[dof,:])
            length(c) == 0 && continue
            push!(related_dofs, c...)
            c = find(G_[dof,:])
            length(c) == 0 && continue
            push!(related_dofs, c...)
        end
        related_dofs = setdiff(related_dofs, all_dofs)
        length(related_dofs) != 0 || return
        info("----- Related equations -----")
        show_rows_in_constraint_matrix(related_dofs, C2, D)
        show_rows_in_constraint_matrix(related_dofs, C2_, D_)
    end

    function print_summary(node_id, dofs)
        s = join(dofs, ", ")
        info()
        info("SUMMARY for node id $node_id with dofs $s:")
        info("----- Current constraint -----")
        info("lambda coefficients in C1 matrix are:")
        show_lambda_coefficients(dofs, C1_)
        info("rows in constraint matrix C2 & D")
        show_rows_in_constraint_matrix(dofs, C2_, D_)
        info("----- New constraint -----")
        info("lambda coefficients in C1 matrix are:")
        show_lambda_coefficients(dofs, C1)
        info("rows in constraint matrix C2 & D")
        show_rows_in_constraint_matrix(dofs, C2, D)
        show_related_equations(dofs, C2, C2_, D, D_)
    end

    info("System is overconstrained by $(length(all_dofs)) dofs.")
    info("Overconstrained_nodes: $(join(nodes, ", ")).")
    info("Following dofs already constrained: $(join(all_dofs, ", ")).")
    for node_id in nodes
        dofs = [2*(node_id-1)+1, 2*(node_id-1)+2]

        show_info && print_summary(node_id, dofs)

        # try to resolve issue automatically
        resolved = false
        for (i, action) in enumerate(actions)
            dofs, resolved = action(node_id, dofs)
            info("algorithm $i solved issue? $resolved")
            if resolved
                break
            end
        end

        if resolved
            info("fixed: new setting is")
            show_info && show_rows_in_constraint_matrix(dofs, C2, D; show_status=false)
            show_info && show_rows_in_constraint_matrix(dofs, C2_, D_; show_status=false)
            show_info && show_related_equations(dofs, C2, C2_, D, D_)
            show_info && info()
            continue
        end

        info("unable to resolve overconstrained situation, not continuing")
        show_rows_in_constraint_matrix(dofs, C2, D; show_status=false)
        show_rows_in_constraint_matrix(dofs, C2_, D_; show_status=false)
        show_related_equations(dofs, C2, C2_, D, D_)
        throw("failed to resolve overconstraint situation")
    end
end

ext: jl


text:
cd(raw"/Users/wapisani/Documents/Programming/Julia/Advent_of_Code/")
open("AoC_2015_6.input") do f
    global lighting_instructions = split(read(f,String),"\n")
end

"""
Returns a tuple of (instruction,start coords, stop coords) from an instruction string
"""
function instruction_parser(instruction::String)
    split_instr = split(instruction," ")
    if split_instr[1] == "toggle"
        instr = "toggle"
        start_coords = split(split_instr[2],",")
        startx = parse(Int,start_coords[1])
        starty = parse(Int,start_coords[2])
        stop_coords = split(split_instr[end],",")
        stopx = parse(Int,stop_coords[1])
        stopy = parse(Int,stop_coords[2])
    else
        instr = split_instr[2]
        start_coords = split(split_instr[3],",")
        startx = parse(Int,start_coords[1])
        starty = parse(Int,start_coords[2])
        stop_coords = split(split_instr[end],",")
        stopx = parse(Int,stop_coords[1])
        stopy = parse(Int,stop_coords[2])
    end

    return (instr,startx,starty,stopx,stopy)
end

global light_grid = zeros(1000,1000)
for instruction in lighting_instructions
    instruction = String(instruction)

    instr_tuple = instruction_parser(instruction)
    light_instr = instr_tuple[1]
    startx = instr_tuple[2]
    starty = instr_tuple[3]
    stopx = instr_tuple[4]
    stopy = instr_tuple[5]

    for i=startx:stopx
        for j=starty:stopy
            # Part 2
            if light_instr == "on"
                light_grid[i,j] += 1
            elseif light_instr == "off"
                if light_grid[i,j] == 0
                    continue
                else
                    light_grid[i,j] -= 1
                end
            elseif light_instr == "toggle"
                light_grid[i,j] += 2
            end
            # Part 1
            # if light_instr == "on"
            #     light_grid[i,j] = 1
            # elseif light_instr == "off"
            #     light_grid[i,j] = 0
            # elseif light_instr == "toggle"
            #     if light_grid[i,j] == 1
            #         light_grid[i,j] = 0
            #     else
            #         light_grid[i,j] = 1
            #     end
            # end
        end
    end
end

#println("There are $(Int(sum(light_grid))) lit lights")
println("The total brightness is $(Int(sum(light_grid)))")


text:
###
  Name     : c5_27_CAPM.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 5/30/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com

 Note #1: go to http://finance.yahoo.com to dowload the last 5-year monthly
          data for IBM and S&P500. The ticker for S&P500 is ^GSPC
      #2: delete the first observation which contains 'null'
      #3: since readtable is deprecated, see c5_31_CAPM.jl which
          uses the CSV.read() function
###

using DataFrames;
using GLM, StatsModels
#
x = readtable("c:/temp/ibmMonthly5years.csv")
p=x[:Adj_Close]
n=length(p)
stockRet=p[2:n]./p[1:(n-1)]-1
#
y = readtable("c:/temp/sp500Monthly5years.csv")
p2=y[:Adj_Close]
n2=length(p2)
mktRet=p2[2:n2]./p2[1:(n2-1)]-1
#
n3=min(length(stockRet),length(mktRet))
data = DataFrame(X=mktRet[1:n3], Y=stockRet[1:n3])
OLS = glm(@formula(Y ~ X), data, Normal(), IdentityLink())


text:
ωcos(N::Integer, i::Integer) = cospi(-2i//N)
ωsin(N::Integer, i::Integer) = sinpi(-2i//N)
ωcos(N, i) = cospi(-2i/N)
ωsin(N, i) = sinpi(-2i/N)

ω(N, i) = ωcos(N, i) + im*ωsin(N, i)

struct fft_meta1
    factors::Vector{Int}
    status::Vector{Int}
    steps::Vector{Int}
    remaining::Vector{Int}
    cumprod::Vector{Int}
    inds::Vector{Int}
end
function fft_meta1(N::Integer)
    factors = factor(Vector, N)
    reverse!(factors)
    status = fill(1,4)#status order is last size, current size, cumulative, it_num
    steps = fill(1, length(factors))
    for i ∈ 2:length(factors) #cumprod, but with a 1 in front.
        steps[i] = factors[i-1] * steps[i-1]
    end
    remaining = N ./ steps
    cp = cumprod(factors)
    fft_meta1(factors, status, steps, remaining, cp, fill(1,length(factors)))

end

function shrink_tree!(x::fft_meta1)
    x.status[3] *= x.status[2]
    x.status[1] = x.status[2]
    x.status[2] = pop!(x.factors)
end

### The plan here is to traverse the remaining tree.
function initial_iterate(x::fft_meta1)
    shrink_tree!(x)
 #   if length(iter) > 0
    iter = x.inds
    resize!(iter, length(x.factors))
    fill!(iter, 1)
    iter[end] = 0
    iter
end
function Base.iterate(x::fft_meta1, iter = initial_iterate(x))
    done(x, iter) && return nothing
    for j ∈ length(iter):-1:1
        if iter[j] == x.factors[j]
            iter[j] = 1
        else
            iter[j] += 1
            break
        end
    end
    iter, iter
end
function done(x::fft_meta1, iter)
    done = true
    for j ∈ eachindex(iter)
        if iter[j] != x.factors[j]
            done = false
            break
        end
    end
    done
end
Base.eltype(::fft_meta1) = Vector{Int}
function Base.length(x::fft_meta1)
    out = 1
    @inbounds for i ∈ 1:length(x.factors)
        out *= x.factors[i]
    end
    out
end

function offset_gap(x::fft_meta1, iter)
    offset = iter[1]
    for i ∈ 2:length(x.factors)
        offset += (iter[i]-1)*x.cumprod[i-1]
    end
    offset, x.cumprod[length(x.factors)]
end

function initial_expr(x::fft_meta1, i, iter, o = :o_, input = :x_)
    expr = quote end
    initial_expr!(expr, x, i, iter, o, input)
end

function initial_expr!(expr, x::fft_meta1, i, iter, output = :o_, input = :x_, ::Type{T} = Float64) where T
#    @show iter
    off, gap = offset_gap(x, iter)
    N = x.status[2]
    initial_expr!(expr, off, gap, N, i, output, input, T)
end
function initial_expr!(expr, off, gap, N, i, output = :o_, input = :x_, ::Type{T} = Float64) where T
    for j ∈ 0:N-1
        push!(expr.args,
                :( $(Symbol( output, 1+i*N+j )) = $(Symbol(input, off)) +
                    $(ω(N, j)) * $(Symbol(input, off+gap)) ) )
        for k ∈ 2:N-1 ##push to new line with += to avoid allocation.
            push!(expr.args,
            :( $(Symbol( output, 1+i*N+j )) +=   $(ω(N, j*k)) * $(Symbol(input, off+gap*k) ) ) )
        end
    end
    expr
end
function gen_initial_expr(fm::fft_meta1, ::Type{T} = Float64) where T
    expr = quote end
    for (i,iter) ∈ enumerate(fm)
        initial_expr!(expr, fm, i-1, iter, :o_, :x_, T)
    end
    expr
end

function next_iter!(expr, fm::fft_meta1, ::Type{T} = Float64) where T
    fm.status[end] += 1
    shrink_tree!(fm)
    l, N, c, it_num = fm.status
    for i ∈ 1:length(fm)
        combine!(expr, l, N, c, i-1, (:o_, :u_)[1 + it_num%2], (:u_, :o_)[1 + it_num%2], T )
    end
end

function combine!(expr, last_N, current_N, cumulative_N, i, input, output, ::Type{T} = Float64) where T
    N = cumulative_N * current_N
    offset = i*N
    ind = 0
    for l ∈ 1:current_N, j ∈ 1:cumulative_N
        push!(expr.args, :( $(Symbol(output, 1+ind+offset )) =
                $(Symbol(input, offset+j)) +
                $(ω(N, ind)) * $(Symbol(input, offset+j+cumulative_N)) ) )
        for k ∈ 2:current_N-1 ##push to new line with += to avoid allocation.
            push!(expr.args,
                :( $(Symbol(output, 1+ind+offset )) += ( $(ω(N, ind*k)) ) * $(Symbol(input, offset+j+cumulative_N*k) )))

        end
        ind += 1
    end
    expr
end

function fft_expression(N::Int, ::Type{T}, ::Type{A}) where {T,A}
    fm = fft_meta1(N)
    expr = gen_initial_expr(fm, T)
    n = length(fm.factors)
    for i ∈ 1:n
        next_iter!(expr, fm, T)
    end
    last_out = (:o_,:u_)[2-fm.status[4]%2]
    out = Expr(:call, A)
    for i ∈ 1:N
        push!(out.args, Symbol( last_out, i ))
    end
    expr = quote @fastmath begin $expr end end
    push!(expr.args, out)
    expr
end

max_issues_repo_name: chriselrod/StaticFFTs

max_forks_repo_name: chriselrod/StaticFFTs


text:
using OrdinaryDiffEq, Test, LinearAlgebra

mutable struct SimType{T} <: DEDataVector{T}
  x::Array{T,1}
  f1::T
end

mutable struct SimType2{T} <: DEDataVector{T}
  x::Vector{T}
  y::Vector{T}
  u::Vector{T}
end

@testset "Data Array Tests (DEDataVector)" begin
  f = function (du,u,p,t)
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
  end

  tstop1 = [5.]
  tstop2 = [8.]
  tstop = [5.;8.]

  condition = function (u,t,integrator)
    t in tstop1
  end

  affect! = function (integrator)
    for c in user_cache(integrator)
      c.f1 = 1.5
    end
  end

  save_positions = (true,true)

  cb = DiscreteCallback(condition, affect!; save_positions=save_positions)

  condition2 = function (u,t,integrator)
    t in tstop2
  end

  affect2! = function (integrator)
    for c in user_cache(integrator)
      c.f1 = -1.5
    end
  end

  save_positions = (true,true)

  cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

  cbs = CallbackSet(cb,cb2)

  u0 = SimType{Float64}([10;10], 0.0)
  prob = ODEProblem(f,u0,(0.0,10.0))
  sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)

  sol(1.5:0.5:2.5)

  @test [sol[i].f1 for i in eachindex(sol)] == [fill(0., 9);1.5*ones(5);-1.5*ones(4)]

  A = Matrix(Diagonal([0.3,0.6,0.9]))
  B = transpose([1 2 3])
  C = [1/3 1/3 1/3]

  function mysystem(t,x,dx,p,u)
      ucalc = u(x,p,t)
      x.u = ucalc
      x.y = C*x.x
      dx .= A*x.x + B*x.u
  end

  input = (x,p,t)->(1*one(t)≤t≤2*one(t) ? [one(t)] : [zero(t)])
  prob = DiscreteProblem((dx,x,p,t)->mysystem(t,x,dx,p,input), SimType2(fill(0., 3), fill(0., 1), fill(0., 1)), (0//1,4//1))
  sln = solve(prob, FunctionMap(scale_by_time=false), dt = 1//10)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, FunctionMap(scale_by_time=true), dt = 1//10)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, Euler(), dt = 1//10)

  @test u1 == [sln[idx].u for idx in 1:length(sln)] # Show that discrete is the same
  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, DP5(), dt = 1//10, adaptive=false)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)
end

######################
# DEDataMatrix
mutable struct SimTypeg{T,T2} <: DEDataMatrix{T}
  x::Array{T,2} # two dimensional
  f1::T2
end

@testset "Data Array Tests (DEDataMatrix)" begin

  tstop1 = [10.0]
  tstop2 = [300.]

  function mat_condition(u,t,integrator)
    t in tstop1
  end

  function mat_condition2(u,t,integrator)
    t in tstop2
  end

  function mat_affect!(integrator)
    for c in user_cache(integrator)
      c.f1 = +1.0
    end
  #  integrator.u[1,1] = 0.001
  end

  function mat_affect2!(integrator)
    for c in user_cache(integrator)
      c.f1 = 0.0
    end
  end

  save_positions = (true,true)
  cb = DiscreteCallback(mat_condition, mat_affect!, save_positions=save_positions)
  save_positions = (false,true)
  cb2 = DiscreteCallback(mat_condition2, mat_affect2!, save_positions=save_positions)
  cbs = CallbackSet(cb,cb2)

  function sigmoid(du,u,p,t)
    du[1,1] = 0.01*u[1,1]*(1-u[1,1]/20)
    du[1,2] = 0.01*u[1,2]*(1-u[1,2]/20)
    du[2,1] = 0.01*u[2,1]*(1-u[2,1]/20)
    du[2,2] = u.f1*du[1,1]
  end

  u0 = SimTypeg(fill(0.00001,2,2),0.0)
  tspan = (0.0,3000.0)
  prob = ODEProblem(sigmoid,u0,tspan)

  tstop =[tstop1;tstop2]
  sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)
  sol = solve(prob,Rodas4(),callback = cbs, tstops=tstop)
  sol = solve(prob,Kvaerno3(),callback = cbs, tstops=tstop)
  @test_broken sol = solve(prob,Rodas4(autodiff=false),callback = cbs, tstops=tstop)
  @test_broken sol = solve(prob,Kvaerno3(autodiff=false),callback = cbs, tstops=tstop)
end


text:
import .Plots.plot
import .Plots.plot!

function plot(env::CartPoleEnv)
    s, a, d = env.state, env.action, env.done
    x, xdot, theta, thetadot = s
    l = 2 * env.params.halflength
    xthreshold = env.params.xthreshold
    # set the frame
    p = plot(
        xlims=(-xthreshold, xthreshold),
        ylims=(-.1, l + 0.1),
        legend=false,
        border=:none,
    )
    # plot the cart
    plot!([x - 0.5, x - 0.5, x + 0.5, x + 0.5], [-.05, 0, 0, -.05];
        seriestype=:shape,
    )
    # plot the pole
    plot!([x, x + l * sin(theta)], [0, l * cos(theta)];
        linewidth=3,
    )
    # plot the arrow
    plot!([x + (a == 1) - 0.5, x + 1.4 * (a == 1)-0.7], [ -.025, -.025];
        linewidth=3,
        arrow=true,
        color=2,
    )
    # if done plot pink circle in top right
    if d
        plot!([xthreshold - 0.2], [l];
            marker=:circle,
            markersize=20,
            markerstrokewidth=0.,
            color=:pink,
        )
    end

    p
end

# adapted from https://github.com/JuliaML/Reinforce.jl/blob/master/src/envs/mountain_car.jl
height(xs) = sin(3 * xs) * 0.45 + 0.55
rotate(xs, ys, θ) = xs * cos(θ) - ys * sin(θ), ys * cos(θ) + xs * sin(θ)
translate(xs, ys, t) = xs .+ t[1], ys .+ t[2]

function plot(env::MountainCarEnv)
    s = env.state
    d = env.done

    p = plot(
        xlims=(env.params.min_pos - 0.1, env.params.max_pos + 0.2),
        ylims=(-.1, height(env.params.max_pos) + 0.2),
        legend=false,
        border=:none,
    )
    # plot the terrain
    xs = LinRange(env.params.min_pos, env.params.max_pos, 100)
    ys = height.(xs)
    plot!(xs, ys)

    # plot the car
    x = s[1]
    θ = cos(3 * x)
    carwidth = 0.05
    carheight = carwidth / 2
    clearance = 0.2 * carheight
    xs = [-carwidth / 2, -carwidth / 2, carwidth / 2, carwidth / 2]
    ys = [0, carheight, carheight, 0]
    ys .+= clearance
    xs, ys = rotate(xs, ys, θ)
    xs, ys = translate(xs, ys, [x, height(x)])
    plot!(xs, ys; seriestype=:shape)

    # if done plot pink circle in top right
    if d
        plot!([xthreshold - 0.2], [l];
            marker=:circle,
            markersize=20,
            markerstrokewidth=0.,
            color=:pink,
        )
    end

    p
 end


text: println(Dates.format(now(), "yyyy-mm-dd"))


text:
using NumericExtensions
using Base.Test

# Test view construction and element accessing

# 1D

a = rand(5)
b = rand(5)
v = unsafe_view(a)

@test ndims(v) == 1
@test length(v) == 5
@test size(v) == (5,)
@test size(v, 1) == 5
@test size(v, 2) == 1
@test pointer(v) == pointer(a)
@test pointer(v, 2) == pointer(a, 2)
@test copy(v) == a

for i in 1 : 5
    @test v[i] == a[i]
    v[i] = b[i]
    @test v[i] == b[i]
end

# 2D

a = rand(4, 5)
b = rand(4, 5)

ac = copy(a)
v = unsafe_view(ac)

@test ndims(v) == 2
@test length(v) == 20
@test size(v) == (4, 5)
@test size(v, 1) == 4
@test size(v, 2) == 5
@test size(v, 3) == 1
@test pointer(v) == pointer(ac)
@test pointer(v, 2) == pointer(ac, 2)
@test copy(v) == a

for i in 1 : 20
    @test v[i] == a[i]
    v[i] = b[i]
    @test v[i] == b[i]
end

ac = copy(a)
v = UnsafeMatrixView{Float64}(pointer(ac), 4, 5)

for j in 1 : 5
    for i in 1 : 4
        @test v[i,j] == a[i,j]
        v[i,j] = b[i,j]
        @test v[i,j] == b[i,j]
    end
end

# 3D

a = rand(4, 5, 3)
b = rand(4, 5, 3)

ac = copy(a)
v = unsafe_view(ac)

@test ndims(v) == 3
@test length(v) == 60
@test size(v) == (4, 5, 3)
@test size(v, 1) == 4
@test size(v, 2) == 5
@test size(v, 3) == 3
@test size(v, 4) == 1
@test pointer(v) == pointer(ac)
@test pointer(v, 2) == pointer(ac, 2)
@test copy(v) == a

for i in 1 : 60
    @test v[i] == a[i]
    v[i] = b[i]
    @test v[i] == b[i]
end

ac = copy(a)
v = unsafe_view(ac)

for j in 1 : 15
    for i in 1 : 4
        @test v[i,j] == a[i,j]
        v[i,j] = b[i,j]
        @test v[i,j] == b[i,j]
    end
end

ac = copy(a)
v = unsafe_view(ac)

for k in 1 : 3
    for j in 1 : 5
        for i in 1 : 4
            @test v[i,j,k] == a[i,j,k]
            v[i,j,k] = b[i,j,k]
            @test v[i,j,k] == b[i,j,k]
        end
    end
end

# Test Subview construction

# 1D

a1 = rand(5)
@test copy(unsafe_view(a1, :)) == a1[:]
@test copy(unsafe_view(a1, :, 1)) == a1[:, 1]
@test copy(unsafe_view(a1, :, 1, 1)) == a1[:, 1, 1]
@test copy(unsafe_view(a1, :, :)) == a1[:, :]
@test copy(unsafe_view(a1, :, :, 1)) == a1[:, :, 1]
@test copy(unsafe_view(a1, :, :, :)) == a1[:, :, :]
@test copy(unsafe_view(a1, :, :, 1:1)) == a1[:, :, 1:1]
@test copy(unsafe_view(a1, :, 1:1)) == a1[:, 1:1]
@test copy(unsafe_view(a1, :, 1:1, 1)) == a1[:, 1:1, 1]

@test copy(unsafe_view(a1, 2:4)) == a1[2:4]
@test copy(unsafe_view(a1, 2:4, 1)) == a1[2:4, 1]
@test copy(unsafe_view(a1, 2:4, 1, 1)) == a1[2:4, 1, 1]

# 2D

a2 = rand(5, 6)
@test copy(unsafe_view(a2, :)) == a2[:]
@test copy(unsafe_view(a2, :, 3)) == a2[:, 3]
@test copy(unsafe_view(a2, :, 3, 1)) == a2[:, 3, 1]
@test copy(unsafe_view(a2, :, :)) == a2[:, :]
@test copy(unsafe_view(a2, :, :, 1)) == a2[:, :, 1]
@test copy(unsafe_view(a2, :, :, :)) == a2[:, :, :]
@test copy(unsafe_view(a2, :, :, 1:1)) == a2[:, :, 1:1]
@test copy(unsafe_view(a2, :, 2:5)) == a2[:, 2:5]
@test copy(unsafe_view(a2, :, 2:5, 1)) == a2[:, 2:5, 1]

@test copy(unsafe_view(a2, 2:4)) == a2[2:4]
@test copy(unsafe_view(a2, 2:4, 3)) == a2[2:4, 3]
@test copy(unsafe_view(a2, 2:4, 3, 1)) == a2[2:4, 3, 1]

# 3D

a3 = rand(5, 6, 3)
@test copy(unsafe_view(a3, :)) == a3[:]
@test copy(unsafe_view(a3, :, 3)) == a3[:, 3]
@test copy(unsafe_view(a3, :, 3, 2)) == a3[:, 3, 2]
@test copy(unsafe_view(a3, :, :)) == a3[:, :]
@test copy(unsafe_view(a3, :, :, 2)) == a3[:, :, 2]
@test copy(unsafe_view(a3, :, :, :)) == a3[:, :, :]
@test copy(unsafe_view(a3, :, :, 1:2)) == a3[:, :, 1:2]
@test copy(unsafe_view(a3, :, 2:5)) == a3[:, 2:5]
@test copy(unsafe_view(a3, :, 2:5, 2)) == a3[:, 2:5, 2]

@test copy(unsafe_view(a3, 2:4)) == a3[2:4]
@test copy(unsafe_view(a3, 2:4, 3)) == a3[2:4, 3]
@test copy(unsafe_view(a3, 2:4, 3, 1)) == a3[2:4, 3, 1]

# 4D

a4 = rand(5, 6, 3, 2)
@test copy(unsafe_view(a4, :)) == a4[:]
@test copy(unsafe_view(a4, :, 3)) == a4[:, 3]
@test copy(unsafe_view(a4, :, 3, 2)) == a4[:, 3, 2]
@test copy(unsafe_view(a4, :, :)) == a4[:, :]
@test copy(unsafe_view(a4, :, :, 2)) == a4[:, :, 2]
@test copy(unsafe_view(a4, :, :, :)) == a4[:, :, :]
@test copy(unsafe_view(a4, :, :, 1:2)) == a4[:, :, 1:2]
@test copy(unsafe_view(a4, :, 2:5)) == a4[:, 2:5]
@test copy(unsafe_view(a4, :, 2:5, 2)) == a4[:, 2:5, 2]

@test copy(unsafe_view(a4, 2:4)) == a4[2:4]
@test copy(unsafe_view(a4, 2:4, 3)) == a4[2:4, 3]
@test copy(unsafe_view(a4, 2:4, 3, 1)) == a4[2:4, 3, 1]


text:
# ### semilog x

#md # ![FILE_NAME.png](images/FILE_NAME.png)

plot(log, 0+0.00001, 5, xscale=:log10, title="semilogx")


text:
"""
"""
struct ExtendedFESpace{S<:SingleFieldFESpace} <: SingleFieldFESpace
  space::S
  model::RestrictedDiscreteModel
  partition::PosNegPartition
  function ExtendedFESpace(space::SingleFieldFESpace,model::RestrictedDiscreteModel)
    model_portion = model.model
    @check get_triangulation(model_portion) === get_triangulation(space)
    @notimplementedif ConstraintStyle(space) == Constrained()
    full_cell_to_cell = Geometry.get_cell_to_parent_cell(model)
    partition = PosNegPartition(full_cell_to_cell,num_cells(get_parent_model(model)))
    new{typeof(space)}(space,model,partition)
  end
end

ConstraintStyle(::Type{ExtendedFESpace{S}}) where S = ConstraintStyle(S)

function get_cell_dof_ids(f::ExtendedFESpace)
  nfull, nvoid = Arrays.pos_and_neg_length(f.partition)
  fullcell_to_ids = get_cell_dof_ids(f.space)
  @check length(fullcell_to_ids) == nfull
  T = eltype(fullcell_to_ids)
  voidcell_to_ids = Fill(similar(T,0),nvoid)
  lazy_map( PosNegReindex(fullcell_to_ids,voidcell_to_ids), f.partition)
end

function get_cell_shapefuns(f::ExtendedFESpace)
  nfull, nvoid = Arrays.pos_and_neg_length(f.partition)
  dv = get_cell_shapefuns(f.space)
  data = get_data(dv)
  fullcell_shapefuns = lazy_map(VoidBasisMap(false),data)
  @check length(fullcell_shapefuns) == nfull
  voidcell_shapefuns = Fill(VoidBasis(testitem(data),true),nvoid)
  a = lazy_map( PosNegReindex(fullcell_shapefuns,voidcell_shapefuns),f.partition)
  SingleFieldFEBasis(a,get_triangulation(f),TestBasis(),DomainStyle(dv))
end

function get_cell_dof_basis(f::ExtendedFESpace)
  nfull, nvoid = Arrays.pos_and_neg_length(f.partition)
  s = get_cell_dof_basis(f.space)
  data = get_data(s)
  fullcell_dof_basis = lazy_map(VoidBasisMap(false),data)
  @check length(fullcell_dof_basis) == nfull
  voidcell_dof_basis = Fill(VoidBasis(testitem(data),true),nvoid)
  a = lazy_map( PosNegReindex(fullcell_dof_basis,voidcell_dof_basis),f.partition)
  CellDof(a,get_triangulation(f),DomainStyle(s))
end

function scatter_free_and_dirichlet_values(f::ExtendedFESpace,fv,dv)
  nfull, nvoid = Arrays.pos_and_neg_length(f.partition)
  fullcell_val = scatter_free_and_dirichlet_values(f.space,fv,dv)
  @check length(fullcell_val) == nfull
  T = eltype(fullcell_val)
  voidcell_val = Fill(similar(T,0),nvoid)
  lazy_map( PosNegReindex(fullcell_val,voidcell_val), f.partition)
end

function gather_free_and_dirichlet_values!(fv,dv,f::ExtendedFESpace,cv)
  _cv = lazy_map(Reindex(cv),f.partition.ipos_to_i)
  gather_free_and_dirichlet_values!(fv,dv,f.space,_cv)
end

# Delegated functions

function get_free_dof_ids(f::ExtendedFESpace)
  get_free_dof_ids(f.space)
end

function get_vector_type(f::ExtendedFESpace)
  get_vector_type(f.space)
end

function get_dirichlet_dof_ids(f::ExtendedFESpace)
  get_dirichlet_dof_ids(f.space)
end

function num_dirichlet_tags(f::ExtendedFESpace)
  num_dirichlet_tags(f.space)
end

function get_dirichlet_dof_tag(f::ExtendedFESpace)
  get_dirichlet_dof_tag(f.space)
end

function get_triangulation(f::ExtendedFESpace)
  get_triangulation(get_parent_model(f.model))
end

# Helpers

struct VoidBasisMap <: Map
  isvoid::Bool
end

@inline Arrays.evaluate!(cache,k::VoidBasisMap,b) = VoidBasis(b,k.isvoid)

struct VoidBasis{T,A} <: AbstractVector{T}
  basis::A
  isvoid::Bool
  function VoidBasis(basis::AbstractVector{T},isvoid::Bool) where T
    new{T,typeof(basis)}(basis,isvoid)
  end
end

function Base.size(a::VoidBasis)
  if a.isvoid
    (0,)
  else
    size(a.basis)
  end
end

Base.IndexStyle(::Type{<:VoidBasis}) = IndexLinear()

function Base.getindex(a::VoidBasis,i::Integer)
  if a.isvoid
    @unreachable "Unable to access 0-length array"
  else
    a.basis[i]
  end
end

Arrays.testitem(a::VoidBasis) = testitem(a.basis)

function Fields.return_cache(a::VoidBasis,x::Point)
  cb = return_cache(a.basis,x)
  bx = return_value(a.basis,x)
  r = similar(bx,(0,))
  cb,r
end

function Fields.return_cache(a::VoidBasis,x::Field)
  cb = return_cache(a.basis,x)
  bx = return_value(a.basis,x)
  r = similar(bx,(0,))
  cb,r
end

function Fields.return_cache(a::VoidBasis,x::AbstractVector{<:Point})
  cb = return_cache(a.basis,x)
  bx = return_value(a.basis,x)
  r = similar(bx,(length(x),0))
  cb,r
end

function Fields.return_cache(a::VoidBasis,v::AbstractVector{<:Field})
  cb = return_cache(a.basis,v)
  bx = return_value(a.basis,v)
  r = similar(bx,(0,length(v)))
  cb,r
end

for T in (:Point,:Field,:(AbstractVector{<:Point}),:(AbstractVector{<:Field}))
  @eval begin

    @inline function Fields.evaluate!(cache,a::VoidBasis,x::$T)
      cb, r = cache
      if a.isvoid
        r
      else
        evaluate!(cb,a.basis,x)
      end
    end

  end
end

function Fields.evaluate!(cache,k::Broadcasting{typeof(∇)},a::VoidBasis)
  VoidBasis(k(a.basis),a.isvoid)
end

function Fields.evaluate!(cache,k::Broadcasting{typeof(∇∇)},a::VoidBasis)
  VoidBasis(k(a.basis),a.isvoid)
end


text:
mutable struct TaggedIntrusiveList{T, Tag}
    head::Union{T, Nothing}

    TaggedIntrusiveList{T, Tag}() where {T, Tag} = new(nothing)
end

const IntrusiveList{T} = TaggedIntrusiveList{T, :default}

@inline function checkbounds(list::IntrusiveList)
    isempty(list) && throw(BoundsError(list))
end

Base.eltype(::Type{TaggedIntrusiveList{T, Tag}}) where {T, Tag} = T

@inline Base.isempty(list::TaggedIntrusiveList) = isnothing(list.head)
@inline Base.empty!(list::TaggedIntrusiveList) = list.head = nothing

function Base.length(list::TaggedIntrusiveList{T, Tag}) where {T, Tag}
    c = 0
    for n in list
        c += 1
    end
    c
end

function Base.show(io::IO, list::TaggedIntrusiveList)
    first = true
    print(io, "[")
    for x in list
        !first && print(io, "↔")
        first = false
        print(io, x)
    end
    print(io, "]")
end

@inline function Base.iterate(list::TaggedIntrusiveList{T, Tag}) where {T, Tag}
    isempty(list) && return nothing
    head = list.head::T
    (head, head)
end

@inline function Base.iterate(list::TaggedIntrusiveList{T, Tag}, state::T) where {T, Tag}
    x = getnext(state, Val{Tag}())::T
    x === list.head ? nothing : (x, x)
end

@inline function Base.first(list::TaggedIntrusiveList{T, Tag}) where {T, Tag}
    @boundscheck checkbounds(list)
    list.head::T
end

@inline function Base.last(list::TaggedIntrusiveList{T, Tag}) where {T, Tag}
    @boundscheck checkbounds(list)
    getprev(list.head::T, Val{Tag}())::T
end

function Base.push!(list::TaggedIntrusiveList{T, Tag}, node::T) where {T, Tag}
    if isempty(list)
        setnextprevself!(node, Val{Tag}())
        list.head = node
    else
        circlist_link_before!(list.head::T, node, Val{Tag}())
    end
    list
end

function Base.pop!(list::TaggedIntrusiveList{T, Tag}) where {T, Tag}
    isempty(list) && throw(ArgumentError("list must be non-empty"))
    head = list.head::T
    tail = getprev(head, Val{Tag}())::T
    if circlist_unlink!(tail, Val{Tag}()) === tail
        list.head = nothing
    end
    tail
end

function Base.pushfirst!(list::TaggedIntrusiveList{T, Tag}, node::T) where {T, Tag}
    if isempty(list)
        setnextprevself!(node, Val{Tag}())
    else
        circlist_link_before!(list.head::T, node, Val{Tag}())
    end
    list.head = node
    list
end

function Base.popfirst!(list::TaggedIntrusiveList{T, Tag}) where {T, Tag}
    isempty(list) && throw(ArgumentError("list must be non-empty"))
    head = list.head::T
    new_head = circlist_unlink!(head, Val{Tag}())
    list.head = head === new_head ? nothing : new_head
    head
end

function Base.delete!(list::TaggedIntrusiveList{T, Tag}, node::T) where {T, Tag}
    node_next = circlist_unlink!(node, Val{Tag}())
    if node === list.head
        if node_next === node
            list.head = nothing
        else
            list.head = node_next
        end
    end
    list
end

function deleteafter!(list::TaggedIntrusiveList{T, Tag}, node::T) where {T, Tag}
    delete!(list, getnext(node, Val{Tag}())::T)
end

function insertafter!(list::TaggedIntrusiveList{T, Tag}, after::T, node::T) where {T, Tag}
    next = getnext(after, Val{Tag}())::T
    circlist_link_before!(next, node, Val{Tag}())
    list
end

function Base.append!(list::TaggedIntrusiveList{T, Tag}, list2::TaggedIntrusiveList{T, Tag}) where {T, Tag}
    isempty(list2) && return list
    if isempty(list)
        list.head = list2.head::T
    else
        head1 = list.head::T
        tail1 = getprev(head1, Val{Tag}())::T
        head2 = list2.head::T
        tail2 = getprev(head2, Val{Tag}())::T
        setnext!(tail1, head2, Val{Tag}())
        setprev!(head2, tail1, Val{Tag}())
        setnext!(tail2, head1, Val{Tag}())
        setprev!(head1, tail2, Val{Tag}())
    end
    list2.head = nothing
    list
end


text:
#=
Created on Monday 13 January 2020
Last update: Saturday 15 August 2020

@author: Michiel Stock
michielfmstock@gmail.com

Test functions for optimization
=#

module TestFuns

using LinearAlgebra

function ackley(x; a=20, b=0.2, c=2π)
    d = length(x)
    return -a * exp(-b*sqrt(sum(x.^2)/d)) -
        exp(sum(cos.(c .* x))/d)
end

ackley(x...; kwargs...) = ackley(x; kwargs...)

# Branin function
function branin((x1, x2); a=1, b=5.1/(4pi^2), c=5/pi, r=6, s=10, t=1/8pi)
    return a * (x2 - b * x1^2 + c * x1 - r)^2 + s * (1 - t) * cos(x1) + s
end

branin(x1, x2; kwargs...) = branin((x1, x2); kwargs...)

# booth
booth((x1, x2)) = (x1+2x2-7)^2 + (2x1+x2-5)^2
booth(x1, x2) = booth((x1, x2))

# Rosenbrock
rosenbrock((x1, x2); a=1, b=5) = (a-x1)^2 + b*(x2-x1^2)^2
rosenbrock(x1, x2; kwargs...) = rosenbrock((x1, x2); kwargs...)

rastrigine(x; A=10) = length(x) * A + sum(x.^2 .+ A .* cos.(2pi .* x))
rastrigine(x...; A=10) = rastrigine(x; A=A)

function flower(x; a=1, b=1, c=4)
    return a * norm(x) + b * sin(c * atan(x[2], x[1]))
end

flower(x1, x2; kwargs...) = flower((x1, x2); kwargs...)

export ackley, branin, rosenbrock, rastrigine, flower, booth

# functions for convex optimization

fquadr((x1, x2); γ=10.0) = 0.5(x1^2 + γ * x2^2);
grad_fquadr((x1, x2); γ=10.0) = [x1, γ * x2];
hess_fquadr((x1, x2); γ=10.0) = [1 0; 0 γ];

export fquadr, grad_fquadr, hess_fquadr

import Zygote  # NOT doing this by hand!

fnonquadr((x1, x2)) = log(exp(x1+3x2-0.1) + exp(x1-3x2-0.1)+exp(-x1-0.1));
grad_fnonquadr(x) = Zygote.gradient(fnonquadr, x)[1];
hess_fnonquadr(x) = Zygote.hessian(fnonquadr, x);

export fnonquadr, grad_fnonquadr, hess_fnonquadr

end  # module TestFuns


text:
function send_logs(cur_player::AbstractPlayer, limit_logs::Bool=true)
  cur_round = cur_player.round

  isempty(cur_round.logs) && return true

  cur_round.is_making_call && return false

  cur_round.is_making_call = true

  end_index = length(cur_round.logs)

  limit_logs && ( end_index = min(30, end_index) )

  cur_log_range = 1:end_index

  cur_logs = cur_round.logs[cur_log_range]

  cur_batch_count = cur_round.batch_count + 1

  cur_request = Dict(
    "round" => Dict(
      "password" => cur_round.client_password,
      "batch_count" => cur_batch_count,
      "logs_attributes" => cur_logs
    )
  )

  is_successful_request = false

  cur_response = nothing

  try
    cur_response = HTTP.put(
      "$(rounds_endpoint)/$(cur_round.server_id)",
      gzip_headers,
      transcode(GzipCompressor, JSON.json(cur_request))
    )

    is_successful_request = ( cur_response.status == 200 )
  catch cur_error
    any(
      cur_error_type -> isa(cur_error, cur_error_type),
      api_errors
    ) || rethrow(cur_error)
  end

  if is_successful_request
    cur_round.batch_count = cur_batch_count

    deleteat!(cur_round.logs, cur_log_range)
  else
    tetris_js(
      cur_player.game.scope,
      JSString("""
        \$(".js-offline").removeClass("hidden");
      """)
    )

    cur_player.round.is_keeping_score = false
  end

  cur_round.is_making_call = false

  true
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# The methods below are specialized on the mortar type
# and called from the basic `create_cache` method at the top.
function create_cache(mesh::P4estMesh{2}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype)
  # TODO: Taal performance using different types
  MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)},
                uEltype, 2,
                nvariables(equations) * nnodes(mortar_l2)}
  fstar_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()]
  fstar_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()]
  u_threaded =           MA2d[MA2d(undef) for _ in 1:Threads.nthreads()]

  (; fstar_upper_threaded, fstar_lower_threaded, u_threaded)
end

#     index_to_start_step_2d(index::Symbol, index_range)
#
# Given a symbolic `index` and an `indexrange` (usually `eachnode(dg)`),
# return `index_start, index_step`, i.e., a tuple containing
# - `index_start`, an index value to begin a loop
# - `index_step`,  an index step to update during a loop
# The resulting indices translate surface indices to volume indices.
#
# !!! warning
#     This assumes that loops using the return values are written as
#
#     i_volume_start, i_volume_step = index_to_start_step_2d(symbolic_index_i, index_range)
#     j_volume_start, j_volume_step = index_to_start_step_2d(symbolic_index_j, index_range)
#
#     i_volume, j_volume = i_volume_start, j_volume_start
#     for i_surface in index_range
#       # do stuff with `i_surface` and `(i_volume, j_volume)`
#
#       i_volume += i_volume_step
#       j_volume += j_volume_step
#     end
@inline function index_to_start_step_2d(index::Symbol, index_range)
  index_begin = first(index_range)
  index_end   = last(index_range)

  if index === :begin
    return index_begin, 0
  elseif index === :end
    return index_end, 0
  elseif index === :i_forward
    return index_begin, 1
  else # if index === :i_backward
    return index_end, -1
  end
end

function prolong2interfaces!(cache, u,
                             mesh::P4estMesh{2},
                             equations, surface_integral, dg::DG)
  @unpack interfaces = cache
  index_range = eachnode(dg)

  @threaded for interface in eachinterface(dg, cache)
    # Copy solution data from the primary element on a case-by-case basis
    # to get the correct face and orientation.
    # Note that in the current implementation, the interface will be
    # "aligned at the primary element", i.e., the index of the primary side
    # will always run forwards.
    primary_element = interfaces.element_ids[1, interface]
    primary_indices = interfaces.node_indices[1, interface]

    i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
    j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)

    i_primary = i_primary_start
    j_primary = j_primary_start
    for i in eachnode(dg)
      for v in eachvariable(equations)
        interfaces.u[1, v, i, interface] = u[v, i_primary, j_primary, primary_element]
      end
      i_primary += i_primary_step
      j_primary += j_primary_step
    end

    # Copy solution data from the secondary element on a case-by-case basis
    # to get the correct face and orientation.
    secondary_element = interfaces.element_ids[2, interface]
    secondary_indices = interfaces.node_indices[2, interface]

    i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], index_range)
    j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], index_range)

    i_secondary = i_secondary_start
    j_secondary = j_secondary_start
    for i in eachnode(dg)
      for v in eachvariable(equations)
        interfaces.u[2, v, i, interface] = u[v, i_secondary, j_secondary, secondary_element]
      end
      i_secondary += i_secondary_step
      j_secondary += j_secondary_step
    end
  end

  return nothing
end

function calc_interface_flux!(surface_flux_values,
                              mesh::P4estMesh{2},
                              nonconservative_terms::Val{false},
                              equations, surface_integral, dg::DG, cache)
  @unpack surface_flux = surface_integral
  @unpack u, element_ids, node_indices = cache.interfaces
  @unpack contravariant_vectors = cache.elements
  index_range = eachnode(dg)

  @threaded for interface in eachinterface(dg, cache)
    # Get information on the primary element, compute the surface fluxes,
    # and store them for the primary element
    primary_element  = element_ids[1, interface]
    primary_indices  = node_indices[1, interface]
    primary_direction = indices2direction(primary_indices)

    i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
    j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)

    i_primary = i_primary_start
    j_primary = j_primary_start
    for i in eachnode(dg)
      u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface)

      # Contravariant vectors at interfaces in negative coordinate direction
      # are pointing inwards. This is handled by `get_normal_direction`.
      normal_direction = get_normal_direction(primary_direction, contravariant_vectors,
                                              i_primary, j_primary, primary_element)
      flux_ = surface_flux(u_ll, u_rr, normal_direction, equations)

      for v in eachvariable(equations)
        surface_flux_values[v, i, primary_direction, primary_element] = flux_[v]
      end

      i_primary += i_primary_step
      j_primary += j_primary_step
    end

    # Get information on the secondary element and copy the numerical fluxes
    # from the primary element to the secondary one
    secondary_element = element_ids[2, interface]
    secondary_indices = node_indices[2, interface]
    secondary_direction = indices2direction(secondary_indices)

    # Note that the index of the primary side will always run forward but
    # the secondary index might need to run backwards for flipped sides.
    if :i_backward in secondary_indices
      for i in eachnode(dg)
        for v in eachvariable(equations)
          surface_flux_values[v, end + 1 - i, secondary_direction, secondary_element] =
            -surface_flux_values[v, i, primary_direction, primary_element]
        end
      end
    else
      for i in eachnode(dg)
        for v in eachvariable(equations)
          surface_flux_values[v, i, secondary_direction, secondary_element] =
            -surface_flux_values[v, i, primary_direction, primary_element]
        end
      end
    end
  end

  return nothing
end

function prolong2boundaries!(cache, u,
                             mesh::P4estMesh{2},
                             equations, surface_integral, dg::DG)
  @unpack boundaries = cache
  index_range = eachnode(dg)

  @threaded for boundary in eachboundary(dg, cache)
    # Copy solution data from the element on a case-by-case basis to get
    # the correct face and orientation.
    element       = boundaries.element_ids[boundary]
    node_indices  = boundaries.node_indices[boundary]

    i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
    j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)

    i_node = i_node_start
    j_node = j_node_start
    for i in eachnode(dg)
      for v in eachvariable(equations)
        boundaries.u[v, i, boundary] = u[v, i_node, j_node, element]
      end
      i_node += i_node_step
      j_node += j_node_step
    end
  end

  return nothing
end

function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing,
                             mesh::P4estMesh{2},
                             equations, surface_integral, dg::DG)
  @unpack boundaries = cache
  @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
  @unpack surface_flux = surface_integral
  index_range = eachnode(dg)

  @threaded for local_index in eachindex(boundary_indexing)
    # Use the local index to get the global boundary index from the pre-sorted list
    boundary = boundary_indexing[local_index]

    # Get information on the adjacent element, compute the surface fluxes,
    # and store them
    element       = boundaries.element_ids[boundary]
    node_indices  = boundaries.node_indices[boundary]
    direction     = indices2direction(node_indices)

    i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
    j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)

    i_node = i_node_start
    j_node = j_node_start
    for i in eachnode(dg)
      # Extract solution data from boundary container
      u_inner = get_node_vars(boundaries.u, equations, dg, i, boundary)

      # Outward-pointing normal direction (not normalized)
      normal_direction = get_normal_direction(direction, contravariant_vectors, i_node, j_node, element)

      # Coordinates at boundary node
      x = get_node_coords(node_coordinates, equations, dg, i_node, j_node, element)

      flux_ = boundary_condition(u_inner, normal_direction, x, t, surface_flux, equations)

      # Copy flux to element storage in the correct orientation
      for v in eachvariable(equations)
        surface_flux_values[v, i, direction, element] = flux_[v]
      end

      i_node += i_node_step
      j_node += j_node_step
    end
  end
end

function prolong2mortars!(cache, u,
                          mesh::P4estMesh{2}, equations,
                          mortar_l2::LobattoLegendreMortarL2,
                          surface_integral, dg::DGSEM)
  @unpack element_ids, node_indices = cache.mortars
  index_range = eachnode(dg)

  @threaded for mortar in eachmortar(dg, cache)
    # Copy solution data from the small elements on a case-by-case basis
    # to get the correct face and orientation.
    small_indices = node_indices[1, mortar]

    i_small_start, i_small_step = index_to_start_step_2d(small_indices[1], index_range)
    j_small_start, j_small_step = index_to_start_step_2d(small_indices[2], index_range)

    for position in 1:2
      i_small = i_small_start
      j_small = j_small_start
      element = element_ids[position, mortar]
      for i in eachnode(dg)
        for v in eachvariable(equations)
          cache.mortars.u[1, v, position, i, mortar] = u[v, i_small, j_small, element]
        end
        i_small += i_small_step
        j_small += j_small_step
      end
    end

    # Buffer to copy solution values of the large element in the correct orientation
    # before interpolating
    u_buffer = cache.u_threaded[Threads.threadid()]

    # Copy solution of large element face to buffer in the
    # correct orientation
    large_indices = node_indices[2, mortar]

    i_large_start, i_large_step = index_to_start_step_2d(large_indices[1], index_range)
    j_large_start, j_large_step = index_to_start_step_2d(large_indices[2], index_range)

    i_large = i_large_start
    j_large = j_large_start
    element = element_ids[3, mortar]
    for i in eachnode(dg)
      for v in eachvariable(equations)
        u_buffer[v, i] = u[v, i_large, j_large, element]
      end
      i_large += i_large_step
      j_large += j_large_step
    end

    # Interpolate large element face data from buffer to small face locations
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 1, :, mortar),
                            mortar_l2.forward_lower,
                            u_buffer)
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 2, :, mortar),
                            mortar_l2.forward_upper,
                            u_buffer)
  end

  return nothing
end

function calc_mortar_flux!(surface_flux_values,
                           mesh::P4estMesh{2},
                           nonconservative_terms::Val{false}, equations,
                           mortar_l2::LobattoLegendreMortarL2,
                           surface_integral, dg::DG, cache)
  @unpack u, element_ids, node_indices = cache.mortars
  @unpack fstar_upper_threaded, fstar_lower_threaded = cache
  @unpack contravariant_vectors = cache.elements
  @unpack surface_flux = surface_integral
  index_range = eachnode(dg)

  @threaded for mortar in eachmortar(dg, cache)
    # Choose thread-specific pre-allocated container
    fstar = (fstar_lower_threaded[Threads.threadid()],
             fstar_upper_threaded[Threads.threadid()])

    # Get information on the small elements, compute the surface fluxes,
    # and store them for the small elements
    small_indices = node_indices[1, mortar]
    small_direction = indices2direction(small_indices)

    i_small_start, i_small_step = index_to_start_step_2d(small_indices[1], index_range)
    j_small_start, j_small_step = index_to_start_step_2d(small_indices[2], index_range)

    # Contravariant vectors at interfaces in negative coordinate direction
    # are pointing inwards. This is handled by `get_normal_direction`.
    for position in 1:2
      i_small = i_small_start
      j_small = j_small_start
      element = element_ids[position, mortar]
      for i in eachnode(dg)
        u_ll, u_rr = get_surface_node_vars(u, equations, dg, position, i, mortar)

        normal_direction = get_normal_direction(small_direction, contravariant_vectors,
                                                i_small, j_small, element)

        flux_ = surface_flux(u_ll, u_rr, normal_direction, equations)

        # Copy flux to buffer
        set_node_vars!(fstar[position], flux_, equations, dg, i)

        i_small += i_small_step
        j_small += j_small_step
      end
    end

    # Buffer to interpolate flux values of the large element to before
    # copying in the correct orientation
    u_buffer = cache.u_threaded[Threads.threadid()]

    mortar_fluxes_to_elements!(surface_flux_values,
                               mesh, equations, mortar_l2, dg, cache,
                               mortar, fstar, u_buffer)
  end

  return nothing
end

@inline function mortar_fluxes_to_elements!(surface_flux_values,
                                            mesh::P4estMesh{2}, equations,
                                            mortar_l2::LobattoLegendreMortarL2,
                                            dg::DGSEM, cache, mortar, fstar, u_buffer)
  @unpack element_ids, node_indices = cache.mortars

  # Copy solution small to small
  small_indices   = node_indices[1, mortar]
  small_direction = indices2direction(small_indices)

  for position in 1:2
    element = element_ids[position, mortar]
    for i in eachnode(dg)
      for v in eachvariable(equations)
        surface_flux_values[v, i, small_direction, element] = fstar[position][v, i]
      end
    end
  end

  # Project small fluxes to large element.
  multiply_dimensionwise!(u_buffer,
                          mortar_l2.reverse_upper, fstar[2],
                          mortar_l2.reverse_lower, fstar[1])

  # The flux is calculated in the outward direction of the small elements,
  # so the sign must be switched to get the flux in outward direction
  # of the large element.
  # The contravariant vectors of the large element (and therefore the normal
  # vectors of the large element as well) are twice as large as the
  # contravariant vectors of the small elements. Therefore, the flux needs
  # to be scaled by a factor of 2 to obtain the flux of the large element.
  u_buffer .*= -2

  # Copy interpolated flux values from buffer to large element face in the
  # correct orientation.
  # Note that the index of the small sides will always run forward but
  # the index of the large side might need to run backwards for flipped sides.
  large_element  = element_ids[3, mortar]
  large_indices  = node_indices[2, mortar]
  large_direction = indices2direction(large_indices)

  if :i_backward in large_indices
    for i in eachnode(dg)
      for v in eachvariable(equations)
        surface_flux_values[v, end + 1 - i, large_direction, large_element] = u_buffer[v, i]
      end
    end
  else
    for i in eachnode(dg)
      for v in eachvariable(equations)
        surface_flux_values[v, i, large_direction, large_element] = u_buffer[v, i]
      end
    end
  end

  return nothing
end

function calc_surface_integral!(du, u,
                                mesh::P4estMesh{2},
                                equations,
                                surface_integral::SurfaceIntegralWeakForm,
                                dg::DGSEM, cache)
  @unpack boundary_interpolation = dg.basis
  @unpack surface_flux_values = cache.elements

  # Note that all fluxes have been computed with outward-pointing normal vectors
  @threaded for element in eachelement(dg, cache)
    for l in eachnode(dg)
      for v in eachvariable(equations)
        # surface at -x
        du[v, 1,          l, element] += (surface_flux_values[v, l, 1, element]
                                          * boundary_interpolation[1, 1])
        # surface at +x
        du[v, nnodes(dg), l, element] += (surface_flux_values[v, l, 2, element]
                                          * boundary_interpolation[nnodes(dg), 2])
        # surface at -y
        du[v, l, 1,          element] += (surface_flux_values[v, l, 3, element]
                                          * boundary_interpolation[1, 1])
        # surface at +y
        du[v, l, nnodes(dg), element] += (surface_flux_values[v, l, 4, element]
                                          * boundary_interpolation[nnodes(dg), 2])
      end
    end
  end

  return nothing
end

end # @muladd


text:
"""
Tools for the exploration of chaos and nonlinear dynamics
"""
module ChaosTools

using DelayEmbeddings
using Entropies
using DynamicalSystemsBase

using DynamicalSystemsBase: DS, DDS, CDS
using DynamicalSystemsBase: MDI, TDI
using DynamicalSystemsBase: stateeltype
using DynamicalSystemsBase.SciMLBase: AbstractODEIntegrator, u_modified!, DEIntegrator
DEI = DEIntegrator

include("deprecations.jl")

include("orbitdiagrams/discrete_diagram.jl")
include("orbitdiagrams/poincare.jl")
include("orbitdiagrams/produce_orbitdiagram.jl")

include("basins/basins_reinit.jl")
include("basins/basins_highlevel.jl")
include("basins/basins_lowlevel.jl")
include("basins/basins_utilities.jl")
include("basins/uncertainty_exp.jl")
include("basins/tipping.jl")

include("dimensions/linear_regions.jl")
include("dimensions/generalized_dim.jl")
include("dimensions/correlationdim.jl")
include("dimensions/correlation_fixedmass.jl")
include("dimensions/molteno.jl")
include("dimensions/kaplanyorke.jl")
include("dimensions/takens_best_estimate.jl")

include("dimreduction/broomhead_king.jl")
include("dimreduction/dyca.jl")

include("period_return/periodic_points.jl")
include("period_return/period.jl")
include("period_return/transit_times_statistics.jl")
include("period_return/fixedpoints.jl")
include("period_return/yin.jl")

include("chaosdetection/lyapunovs.jl")
include("chaosdetection/gali.jl")
include("chaosdetection/expansionentropy.jl")
include("chaosdetection/partially_predictable.jl")
include("chaosdetection/testchaos01.jl")

end # module