text:
module Selection

using ..Ahorn, Maple

displayName = "Selection"
group = "Placements"

drawingLayers = Ahorn.Layer[]

toolsLayer = nothing
targetLayer = nothing

# Drag selection, track individually for easy replacement
selectionRect = Ahorn.Rectangle(0, 0, 0, 0)
selectionPreviews = Set{Ahorn.SelectedObject}()
selections = Set{Ahorn.SelectedObject}()

lastX, lastY = -1, -1
shouldDrag = false
canDrag = false

decalScaleVals = (1.0, 2.0^4)

canResize = false
resizeCursor = "default"
cursorResizeThreshold = 1

const resizeCursorDirections = Dict{String, Tuple{Int, Int}}(
    "n-resize" => (0, -1),
    "ne-resize" => (1, -1),
    "e-resize" => (1, 0),
    "se-resize" => (1, 1),
    "s-resize" => (0, 1),
    "sw-resize" => (-1, 1),
    "w-resize" => (-1, 0),
    "nw-resize" => (-1, -1),
)

relevantRoom = Ahorn.loadedState.room
selectionsClipboard = Set{Ahorn.SelectedObject}()

function copySelections(cut::Bool=false)
    if isempty(selections)
        return false
    end

    empty!(selectionsClipboard)
    union!(selectionsClipboard, deepcopy(selections))

    if cut
        handleDeletion(selections)
        empty!(selections)
    end

    toolsLayer.redraw = true
    redrawTargetLayer!(targetLayer, selectionsClipboard)

    return cut
end

cutSelections() = copySelections(true)

function getCursorForResize(x::Int, y::Int, rect::Ahorn.Rectangle, canHorizontal::Bool=true, canVertical::Bool=true, threshold::Number=cursorResizeThreshold)
    if !(canHorizontal || canVertical)
        return "default"
    end

    res = ""

    onHorizontal = rect.x - threshold <= x <= rect.x + rect.w + threshold
    onVertical = rect.y - threshold <= y <= rect.y + rect.h + threshold

    onNorthLine = onHorizontal && (rect.y - threshold <= y <= rect.y + threshold)
    onSouthLine = onHorizontal && (rect.y + rect.h - threshold <= y <= rect.y + rect.h + threshold)
    onWestLine = onVertical && (rect.x - threshold <= x <= rect.x + threshold)
    onEastLine = onVertical && (rect.x + rect.w - threshold <= x <= rect.x + rect.w + threshold)

    if onNorthLine && canVertical
        res *= "n"

    elseif onSouthLine && canVertical
        res *= "s"
    end

    if onWestLine && canHorizontal
        res *= "w"

    elseif onEastLine && canHorizontal
        res *= "e"
    end

    return isempty(res) ? "default" : res * "-resize"
end

# Windows doesn't have "grabbing" and "grab"
function updateCursor()
    cursor = "default"

    if canResize
        cursor = resizeCursor

    elseif shouldDrag
        cursor = Sys.iswindows() ? "default" : "grabbing"

    elseif canDrag
        cursor = Sys.iswindows() ? "default" : "grab"
    end

    Ahorn.setWindowCursor!(Ahorn.window, cursor)
end

function getSelectionLeftCorner(selections::Set{Ahorn.SelectedObject})
    tlx, tly = typemax(Int), typemax(Int)

    for selection in selections
        x, y = Ahorn.position(selection.target)

        tlx = min(tlx, x)
        tly = min(tly, y)
    end

    return tlx, tly
end

function pasteSelections()
    if isempty(selectionsClipboard)
       return false
    end

    newSelections = deepcopy(selectionsClipboard)
    room = relevantRoom

    offsetX, offsetY = 0, 0

    if Ahorn.cursor !== nothing
        topLeftX, topLeftY = getSelectionLeftCorner(newSelections)

        offsetX = -topLeftX + Ahorn.cursor.x * 8 - 8
        offsetY = -topLeftY + Ahorn.cursor.y * 8 - 8
    end

    for selection in newSelections
        layer, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node

        if layer == "fgDecals"
            target.x += offsetX
            target.y += offsetY

            push!(room.fgDecals, target)

        elseif layer == "bgDecals"
            target.x += offsetX
            target.y += offsetY

            push!(room.bgDecals, target)

        elseif layer == "entities" && node == 0
            target.x += offsetX
            target.y += offsetY

            if haskey(target.data, "nodes")
                target.nodes = Tuple{Integer, Integer}[(x + offsetX, y + offsetY) for (x, y) in target.nodes]
            end

            target.id = Ahorn.EntityIds.nextId()

            push!(room.entities, target)

        elseif layer == "triggers" && node == 0
            target.x += offsetX
            target.y += offsetY

            if haskey(target.data, "nodes")
                target.nodes = Tuple{Integer, Integer}[(x + offsetX, y + offsetY) for (x, y) in target.nodes]
            end

            target.id = Ahorn.EntityIds.nextId()

            push!(room.triggers, target)

        elseif layer == "fgTiles" || layer == "bgTiles"
            target.offsetX += offsetX
            target.offsetY += offsetY

            target.selection = Ahorn.Rectangle(
                target.selection.x + offsetX,
                target.selection.y + offsetY,
                target.selection.w,
                target.selection.h
            )
        end
    end

    finalizeSelections!(selections)
    empty!(selections)
    union!(selections, newSelections)

    toolsLayer.redraw = true
    redrawTargetLayer!(targetLayer, selections)

    return true
end

hotkeys = Ahorn.Hotkey[
    Ahorn.Hotkey(
        "ctrl + c",
        copySelections
    ),
    Ahorn.Hotkey(
        "ctrl + x",
        cutSelections
    ),
    Ahorn.Hotkey(
        "ctrl + v",
        pasteSelections
    )
]

function drawSelections(layer::Ahorn.Layer, room::Ahorn.DrawableRoom, camera::Ahorn.Camera)
    drawnTargets = Set()
    ctx = Ahorn.getSurfaceContext(toolsLayer.surface)

    # Prone to errors with weird pointers, bandaid fix for now
    # This might not actually be needed, but i'm keeping it for now
    if ctx.surface.ptr == C_NULL || ctx.ptr == C_NULL
        return true
    end

    if selectionRect !== nothing && selectionRect.w > 0 && selectionRect.h > 0 && !shouldDrag
        Ahorn.drawRectangle(ctx, selectionRect, Ahorn.colors.selection_selection_fc, Ahorn.colors.selection_selection_bc)
    end

    # Make sure fgTiles render after the bgTiles
    # Looks better, shouldn't cost to much performance
    selectionsArray = collect(selections)
    sort!(selectionsArray, by=r -> r.layerName)

    for selection in selectionsArray
        layer, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node

        if isa(target, Maple.Entity) && !(target in [row[1] for row in drawnTargets])
            Ahorn.renderEntitySelection(ctx, toolsLayer, target, relevantRoom)
        end

        if isa(target, Maple.Trigger) && !(target in [row[1] for row in drawnTargets])
            Ahorn.renderTriggerSelection(ctx, toolsLayer, target, relevantRoom)
        end

        if isa(target, Ahorn.TileSelection)
            Ahorn.drawFakeTiles(ctx, relevantRoom, target.tiles, relevantRoom.objTiles, target.fg, target.selection.x, target.selection.y, clipEdges=true)
        end

        push!(drawnTargets, (target, node))

        # Get a new selection rectangle
        # This is easier than editing the existing rect
        rect = Ahorn.getSelection(target, room.room)
        if isa(rect, Array{Ahorn.Rectangle}) && length(rect) >= node + 1
            Ahorn.drawRectangle(ctx, rect[node + 1], Ahorn.colors.selection_selected_fc, Ahorn.colors.selection_selected_bc)

        else
            Ahorn.drawRectangle(ctx, rect, Ahorn.colors.selection_selected_fc, Ahorn.colors.selection_selected_bc)
        end
    end

    for preview in selectionPreviews
        layer, box, target, node = preview.layerName, preview.rectangle, preview.target, preview.node

        if isa(target, Maple.Entity) && !(target in [row[1] for row in drawnTargets])
            Ahorn.renderEntitySelection(ctx, toolsLayer, target, relevantRoom)
        end

        if isa(target, Maple.Trigger) && !(target in [row[1] for row in drawnTargets])
            Ahorn.renderTriggerSelection(ctx, toolsLayer, target, relevantRoom)
        end

        if !((target, node) in drawnTargets) && !(preview in selections)
            Ahorn.drawRectangle(ctx, box, Ahorn.colors.selection_preview_fc, Ahorn.colors.selection_preview_bc)
        end

        push!(drawnTargets, (target, node))
    end

    return true
end

function clearDragging!()
    global lastX = -1
    global lastY = -1
    global shouldDrag = false
    global canDrag = false
end

function clearResize!()
    global resizeCursor = "default"
    global canResize = false
end

function cleanup()
    finalizeSelections!(selections)
    empty!(selections)

    global selectionRect = Ahorn.Rectangle(0, 0, 0, 0)
    global relevantRoom = Ahorn.loadedState.room
    global targetLayer = nothing
    global toolsLayer = nothing
    global drawingLayers = Ahorn.Layer[]

    clearDragging!()
    clearResize!()
    updateCursor()
end

function toolSelected(subTools::Ahorn.ListContainer, layers::Ahorn.ListContainer, materials::Ahorn.ListContainer)
    global relevantRoom = Ahorn.loadedState.room

    wantedLayer = get(Ahorn.persistence, "placements_layer", "entities")
    Ahorn.updateLayerList!(vcat(["all"], Ahorn.selectableLayers), row -> row[1] == Ahorn.layerName(targetLayer))

    Ahorn.redrawingFuncs["tools"] = drawSelections
    Ahorn.redrawLayer!(toolsLayer)
end

function layerSelected(list::Ahorn.ListContainer, materials::Ahorn.ListContainer, selected::String)
    global targetLayer = Ahorn.getLayerByName(drawingLayers, selected)
    Ahorn.persistence["placements_layer"] = selected
end

function updateSelectionPreviews(x::Int, y::Int)
    if get(Ahorn.config, "preview_possible_selections", true)
        if selectionRect == Ahorn.Rectangle(0, 0, 0, 0)
            # Don't preview tile selections or "all" selection
            # It looks weird and it can cause deletion of tiles

            name = Ahorn.layerName(targetLayer)
            if name != "fgTiles" && name != "bgTiles" && name != "all"
                previousPreviews = deepcopy(selectionPreviews)
                finalizeSelections!(selectionPreviews)
                empty!(selectionPreviews)
                rect = Ahorn.Rectangle(x, y, 1, 1)
                properlyUpdateSelections!(rect, selectionPreviews, best=true)

                # TODO - Unhaunt haunted branch, this works for now
                #sameSelections = previousPreviews == selectionPreviews
                sameSelections = length(previousPreviews) == length(selectionPreviews) && sort(collect(previousPreviews)) == sort(collect(selectionPreviews))

                if !sameSelections
                    Ahorn.redrawLayer!(toolsLayer)
                end
            end
        end
    end
end

function updateResize(x::Int, y::Int)
    if !shouldDrag && length(selections) == 1 && !Ahorn.mouseButtonHeld(0x1)
        selection = first(selections)
        layer, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node
        selectionRect = Ahorn.getSelection(target, Ahorn.loadedState.room)
        rect = isa(selectionRect, Array{Ahorn.Rectangle, 1}) ? selectionRect[1] : selectionRect

        if isa(rect, Ahorn.Rectangle) && (isa(target, Maple.Entity) || isa(target, Maple.Trigger)) && node == 0
            resizeW, resizeH = Ahorn.canResizeWrapper(target)

            name = getCursorForResize(x, y, rect, resizeW, resizeH)

            global resizeCursor = name
            global canResize = name != "default"
        end
    end
end

function updateDrag(x::Int, y::Int)
    if !shouldDrag
        target = Ahorn.hasSelectionAt(selections, Ahorn.Rectangle(x, y, 1, 1), relevantRoom)
        global canDrag = target != false
    end
end

function selectionMotionAbs(rect::Ahorn.Rectangle)
    if rect != selectionRect && !shouldDrag && !canResize
        if get(Ahorn.config, "preview_possible_selections", true)
            # Don't preview tile selections or "all" selection
            # It looks weird and it can cause deletion of tiles

            name = Ahorn.layerName(targetLayer)
            if name != "fgTiles" && name != "bgTiles" && name != "all"
                finalizeSelections!(selectionPreviews)
                empty!(selectionPreviews)

                properlyUpdateSelections!(rect, selectionPreviews)
            end
        end

        global selectionRect = rect

        Ahorn.redrawLayer!(toolsLayer)
    end
end

function selectionMotionAbs(x1::Number, y1::Number, x2::Number, y2::Number)
    ctrl = Ahorn.modifierControl()

    if lastX == -1 && lastY == -1
        global lastX = ctrl ? x1 : div(x1, 8) * 8
        global lastY = ctrl ? y1 : div(y1, 8) * 8

        if !shouldDrag && canDrag
            global shouldDrag = true

            Ahorn.History.addSnapshot!(Ahorn.History.MultiSnapshot("Selections", Ahorn.History.Snapshot[
                Ahorn.History.RoomSnapshot("Selections", Ahorn.loadedState.room),
                Ahorn.History.SelectionSnapshot("Selections", relevantRoom, selections)
            ]))
        end
    end

    if !Ahorn.modifierControl()
        x1 = div(x1, 8) * 8
        y1 = div(y1, 8) * 8

        x2 = div(x2, 8) * 8
        y2 = div(y2, 8) * 8
    end

    dx = x2 - lastX
    dy = y2 - lastY

    global lastX = x2
    global lastY = y2

    if canResize
        if (dx != 0 || dy != 0) && length(selections) == 1
            selection = first(selections)
            layer, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node

            if isa(target, Maple.Trigger) || isa(target, Maple.Entity)
                changed = false
                resizeX, resizeY = get(resizeCursorDirections, resizeCursor, (0, 0))

                resizeWidth, resizeHeight = Ahorn.canResizeWrapper(target)
                minimumWidth, minimumHeight = Ahorn.minimumSizeWrapper(target)

                if resizeX != 0 && resizeWidth
                    if target.width + dx * resizeX >= minimumWidth
                        target.width = target.width + dx * resizeX
                        target.x = target.x + dx * (resizeX < 0)

                        changed |= true
                    end
                end

                if resizeY != 0 && resizeHeight
                    if target.height + dy * resizeY >= minimumHeight
                        target.height = target.height + dy * resizeY
                        target.y = target.y + dy * (resizeY < 0)

                        changed |= true
                    end
                end

                if changed
                    snapshot = Ahorn.History.MultiSnapshot("Selections", Ahorn.History.Snapshot[
                        Ahorn.History.RoomSnapshot("Selections", Ahorn.loadedState.room),
                        Ahorn.History.SelectionSnapshot("Selections", relevantRoom, selections)
                    ])

                    Ahorn.History.addSnapshot!(snapshot)

                    toolsLayer.redraw = true
                    Ahorn.redrawLayer!(Ahorn.getLayerByName(drawingLayers, layer))
                end
            end
        end

    elseif shouldDrag
        global lastX = ctrl ? x2 : div(x2, 8) * 8
        global lastY = ctrl ? y2 : div(y2, 8) * 8

        if dx != 0 || dy != 0
            for selection in selections
               layer, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node

                if applyMovement!(target, dx, dy, node)
                    notifyMovement!(target)
                end
            end

            toolsLayer.redraw = true
            redrawTargetLayer!(targetLayer, selections, String["fgTiles", "bgTiles"], onlyMark=true)
            Ahorn.redrawCanvas!()
        end
    end
end

function mouseMotionAbs(x::Number, y::Number)
    updateDrag(x, y)
    updateSelectionPreviews(x, y)
    updateResize(x, y)

    updateCursor()
end

function properlyUpdateSelections!(rect::Ahorn.Rectangle, selections::Set{Ahorn.SelectedObject}; best::Bool=false, mass::Bool=false)
    retain = Ahorn.modifierShift()

    # Do this before we get new selections
    # This way tiles are settled back into place before we select
    if !retain
        finalizeSelections!(selections)
    end

    unselected, newlySelected = Ahorn.updateSelections!(selections, relevantRoom, Ahorn.layerName(targetLayer), rect, retain=retain, best=best, mass=mass)
    initSelections!(newlySelected)
end

function properlyMassSelection!(selections::Set{Ahorn.SelectedObject}, rect::Ahorn.Rectangle; strict::Bool=false)
    retain = Ahorn.modifierShift()

    # Do this before we get new selections
    # This way tiles are settled back into place before we select
    if !retain
        finalizeSelections!(selections)
    end

    unselected, newlySelected = Ahorn.updateSelections!(selections, relevantRoom, Ahorn.layerName(targetLayer), rect, retain=retain, best=strict, mass=true)
    initSelections!(newlySelected)
end

function getLayersSelected(selections::Set{Ahorn.SelectedObject})
    return unique([selection.layerName for selection in selections])
end

function redrawTargetLayer!(layer::Ahorn.Layer, selections::Set{Ahorn.SelectedObject}, ignore::Array{String, 1}=String[]; onlyMark::Bool=false)
    redrawTargetLayer!(layer, getLayersSelected(selections), ignore, onlyMark=onlyMark)
end

function redrawTargetLayer!(layer::Ahorn.Layer, layers::Array{String, 1}, ignore::Array{String, 1}=String[]; onlyMark::Bool=false)
    needsRedraw = filter(v -> !(v in ignore), layers)

    for layer in needsRedraw
        Ahorn.getLayerByName(drawingLayers, layer).redraw = true
    end

    if !isempty(needsRedraw) || onlyMark
        Ahorn.redrawCanvas!()
    end
end

function selectionFinishAbs(rect::Ahorn.Rectangle)
    # If we are draging we are techically not making a new selection
    if !shouldDrag
        properlyUpdateSelections!(rect, selections)
    end

    clearDragging!()
    clearResize!()
    updateCursor()

    global selectionRect = Ahorn.Rectangle(0, 0, 0, 0)

    Ahorn.redrawLayer!(toolsLayer)
end

function leftClickAbs(x::Number, y::Number)
    rect = Ahorn.Rectangle(x, y, 1, 1)
    properlyUpdateSelections!(rect, selections, best=true)

    clearDragging!()
    clearResize!()
    updateCursor()

    Ahorn.redrawLayer!(toolsLayer)
end

function doubleLeftClickAbs(x::Number, y::Number)
    strict = Ahorn.modifierControl()
    rect = Ahorn.Rectangle(x, y, 1, 1)
    properlyMassSelection!(selections, rect, strict=strict)

    clearDragging!()
    clearResize!()

    Ahorn.redrawLayer!(toolsLayer)
end

function rightClickAbs(x::Number, y::Number)
    if !shouldDrag
        Ahorn.displayProperties(x, y, relevantRoom, targetLayer, toolsLayer, selections)

        clearDragging!()
        clearResize!()

        Ahorn.redrawLayer!(toolsLayer)
    end
end

function layersChanged(layers::Array{Ahorn.Layer, 1})
    wantedLayer = get(Ahorn.persistence, "placements_layer", "entities")

    global drawingLayers = layers
    global toolsLayer = Ahorn.getLayerByName(layers, "tools")
    global targetLayer = Ahorn.selectLayer!(layers, wantedLayer, "entities")
end

function applyTileSelecitonBrush!(target::Ahorn.TileSelection, clear::Bool=false)
    Maple.updateTileSize!(relevantRoom, '0', '0')

    roomTiles = target.fg ? relevantRoom.fgTiles : relevantRoom.bgTiles
    tiles = clear ? fill('0', size(target.tiles)) : target.tiles

    x, y = floor(Int, target.selection.x / 8), floor(Int, target.selection.y / 8)
    brush = Ahorn.Brush(
        "Selection Finisher",
        clear ? fill(1, size(tiles) .- 2) : tiles[2:end - 1, 2:end - 1] .!= '0'
    )

    Ahorn.applyBrush!(brush, roomTiles, tiles[2:end - 1, 2:end - 1], x + 1, y + 1)
end

function afterUndo(map::Maple.Map)
    global selections = Ahorn.fixSelections(relevantRoom, selections)
    Ahorn.redrawLayer!(toolsLayer)
end

function afterRedo(map::Maple.Map)
    global selections = Ahorn.fixSelections(relevantRoom, selections)
    Ahorn.redrawLayer!(toolsLayer)
end

function finalizeSelections!(targets::Set{Ahorn.SelectedObject})
    shouldRedraw = false

    for selection in targets
        layer = selection.layerName

        if layer == "fgTiles" || layer == "bgTiles"
            applyTileSelecitonBrush!(selection.target, false)
            shouldRedraw = true
        end
    end

    if shouldRedraw
        redrawTargetLayer!(targetLayer, targets)
    end
end

function initSelections!(targets::Set{Ahorn.SelectedObject})
    shouldRedraw = false

    for selection in targets
        layer = selection.layerName

        if layer == "fgTiles" || layer == "bgTiles"
            applyTileSelecitonBrush!(selection.target, true)
            shouldRedraw = true
        end
    end

    if shouldRedraw
        redrawTargetLayer!(targetLayer, targets)
    end
end

function setSelections(map::Maple.Map, room::Maple.Room, newSelections::Set{Ahorn.SelectedObject})
    if room.name == relevantRoom.name
        empty!(selections)
        union!(selections, Ahorn.fixSelections(relevantRoom, newSelections))

        Ahorn.redrawLayer!(toolsLayer)
    end
end

function getSelections()
    return selections
end

function gridSnapped(v, step, direction)
    if direction == 0
        return v

    elseif direction > 0
        return v + step - mod(v, step)

    else
        return v - mod1(v, step)
    end
end

applyMovement!(target, ox, oy, node) = false

function applyMovement!(target::Union{Maple.Entity, Maple.Trigger}, ox, oy, node=0)
    if node == 0
        target.data["x"] += ox
        target.data["y"] += oy

    else
        nodes = get(target.data, "nodes", ())

        if length(nodes) >= node
            nodes[node] = nodes[node] .+ (ox, oy)
        end
    end

    return true
end

function applyMovement!(decal::Maple.Decal, ox, oy, node=0)
    decal.x += ox
    decal.y += oy

    return true
end

function applyMovement!(target::Ahorn.TileSelection, ox, oy, node=0)
    target.offsetX += ox
    target.offsetY += oy

    target.selection = Ahorn.Rectangle(target.startX + floor(target.offsetX / 8) * 8, target.startY + floor(target.offsetY / 8) * 8, target.selection.w, target.selection.h)

    return true
end

function applyGridMovement!(target::Union{Maple.Entity, Maple.Trigger}, gridSize, directionX, directionY, node=0)
    if node == 0
        x, y = target.data["x"], target.data["y"]
        ox, oy = gridSnapped(x, gridSize, directionX), gridSnapped(y, gridSize, directionY)

        applyMovement!(target, ox - x, oy - y, node)

    else
        nodes = get(target.data, "nodes", ())

        if length(nodes) >= node
            nx, ny = nodes[node]
            ox, oy = gridSnapped(nx, gridSize, directionX), gridSnapped(ny, gridSize, directionY)

            applyMovement!(target, ox - nx, oy - ny, node)
        end
    end
end

function applyGridMovement!(decal::Maple.Decal, gridSize, directionX, directionY, node=0)
    ox, oy = gridSnapped(decal.x, gridSize, directionX), gridSnapped(decal.y, gridSize, directionY)

    applyMovement!(decal, ox - decal.x, oy - decal.y)
end

function applyGridMovement!(target::Ahorn.TileSelection, gridSize, directionX, directionY, node=0)
    ox, oy = gridSnapped(target.offsetX, gridSize, directionX), gridSnapped(target.offsetY, gridSize, directionY)

    applyMovement!(target, ox - target.offsetX, oy - target.offsetY)
end

function notifyMovement!(entity::Maple.Entity)
    Ahorn.eventToModules(Ahorn.loadedEntities, "moved", entity)
    Ahorn.eventToModules(Ahorn.loadedEntities, "moved", entity, relevantRoom)
end

function notifyMovement!(trigger::Maple.Trigger)
    Ahorn.eventToModules(Ahorn.loadedTriggers, "moved", trigger)
    Ahorn.eventToModules(Ahorn.loadedTriggers, "moved", trigger, relevantRoom)
end

function notifyMovement!(decal::Maple.Decal)
    # Decals doesn't care
end

function notifyMovement!(target::Ahorn.TileSelection)
    # Tiles doesn't care
end

resizeModifiers = Dict{Integer, Tuple{Number, Number}}(
    # w, h
    # Decrease / Increase width
    Int('q') => (1, 0),
    Int('e') => (-1, 0),

    # Decrease / Increase height
    Int('a') => (0, 1),
    Int('d') => (0, -1)
)

addNodeKeys = [Int('n'), Int('+')]

# Turns out having scales besides -1 and 1 on decals causes weird behaviour?
scaleMultipliers = Dict{Integer, Tuple{Number, Number}}(
    # Vertical Flip
    Int('v') => (1, -1),

    # Horizontal Flip
    Int('h') => (-1, 1),
)

# Consider exposing the grid snap value
function handleMovement(event::Ahorn.eventKey)
    redraw = false
    step = Ahorn.modifierControl() ? 1 : 8
    snapMode = get(Ahorn.config, "use_grid_snapping", true)

    for selection in selections
        name, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node
        dirX, dirY = Ahorn.moveDirections[event.keyval]

        if snapMode
            if applyGridMovement!(target, step, dirX, dirY, node)
                notifyMovement!(target)

                redraw = true
            end

        else
            redraw = applyMovement!(target, dirX * step, dirY * step, node)

            if redraw
                notifyMovement!(target)
            end
        end
    end

    return redraw
end

function handleResize(event::Ahorn.eventKey)
    redraw = false
    step = Ahorn.modifierControl() ? 1 : 8

    processed = Set{Union{Maple.Entity, Maple.Trigger}}()
    snapMode = get(Ahorn.config, "use_grid_snapping", true)

    for selection in selections
        name, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node
        dirW, dirH = resizeModifiers[event.keyval]

        if (name == "entities" || name == "triggers") && !(target in processed)
            horizontal, vertical = Ahorn.canResizeWrapper(target)
            minWidth, minHeight = Ahorn.minimumSizeWrapper(target)

            baseWidth = get(target.data, "width", minWidth)
            baseHeight = get(target.data, "height", minHeight)

            if snapMode
                target.data["width"] = horizontal ? max(gridSnapped(baseWidth, step, dirW), minWidth) : baseWidth
                target.data["height"] = vertical ? max(gridSnapped(baseHeight, step, dirH), minHeight) : baseHeight

            else
                target.data["width"] = horizontal ? max(baseWidth + dirW * step, minWidth) : baseWidth
                target.data["height"] = vertical ? max(baseHeight + dirH * step, minHeight) : baseHeight
            end

            redraw = true

            push!(processed, target)

        elseif name == "fgDecals" || name == "bgDecals"
            extraW, extraH = resizeModifiers[event.keyval]
            minVal, maxVal = decalScaleVals

            target.scaleX = sign(target.scaleX) * clamp(abs(target.scaleX) * 2.0^extraW, minVal, maxVal)
            target.scaleY = sign(target.scaleY) * clamp(abs(target.scaleY) * 2.0^extraH, minVal, maxVal)

            redraw = true
        end
    end

    return redraw
end

# Requires a runtime dispatch regardless, using isa for less boilerplate
function handleScaling(event::Ahorn.eventKey)
    redraw = false
    for selection in selections
        msx, msy = scaleMultipliers[event.keyval]
        target = selection.target

        if isa(target, Maple.Decal)
            target.scaleX *= msx
            target.scaleY *= msy

            redraw = true
        end
    end

    return redraw
end

function handleAddNodes(event::Ahorn.eventKey)
    redraw = false

    for selection in selections
        name, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node

        if name == "entities" || name == "triggers"
            least, most = Ahorn.nodeLimits(target)
            nodes = get(target.data, "nodes", Tuple{Int, Int}[])

            if most == -1 || length(nodes) + 1 <= most
                x, y = target.data["x"], target.data["y"]

                if node > 0
                    x, y = nodes[node]
                end

                insert!(nodes, node + 1, (x + 16, y))
                redraw = true

                target.data["nodes"] = nodes
            end
        end
    end

    return redraw
end

function handleDeletion(selections::Set{Ahorn.SelectedObject})
    res = !isempty(selections)
    selectionsArray = collect(selections)

    # Split into different arrays
    tileSelections = filter(s -> s.layerName == "bgTiles" || s.layerName == "fgTiles", selectionsArray)
    entityTriggerSelections = filter(s -> s.layerName == "entities" || s.layerName == "triggers", selectionsArray)
    decalSelections = filter(s -> s.layerName == "bgDecals" || s.layerName == "fgDecals", selectionsArray)

    # Sort entities, otherwise deletion will break with nodes
    sort!(entityTriggerSelections, by=r -> (r.target.id, r.node), rev=true)

    # Deletion for entities and triggers
    for selection in entityTriggerSelections
        name, box, target, node = selection.layerName, selection.rectangle, selection.target, selection.node
        targetList = Ahorn.selectionTargets[name](relevantRoom)

        index = findfirst(isequal(target), targetList)
        if index !== nothing
            if node == 0
                deleteat!(targetList, index)

            else
                least, most = Ahorn.nodeLimits(target)
                nodes = get(target.data, "nodes", [])

                # Delete the node if that doesn't result in too few nodes
                # Delete the whole entity if it does
                if length(nodes) - 1 >= least && length(nodes) >= node
                    deleteat!(nodes, node)

                else
                    deleteat!(targetList, index)
                end
            end
        end
    end

    # Deletion for decals
    for selection in decalSelections
        targetList = Ahorn.selectionTargets[selection.layerName](relevantRoom)

        index = findfirst(isequal(selection.target), targetList)
        if index !== nothing
            deleteat!(targetList, index)
        end
    end

    # Tiles are deleted by removing them from the set, no special handle

    if !isempty(selections)
        empty!(selections)
    end

    return res
end

# Refactor and prettify code once we know how to handle tiles here,
# this also includes the handle functions
function keyboard(event::Ahorn.eventKey)
    needsRedraw = false
    layersSelected = getLayersSelected(selections)
    snapshot = Ahorn.History.MultiSnapshot("Selections", Ahorn.History.Snapshot[
        Ahorn.History.RoomSnapshot("Selections", Ahorn.loadedState.room),
        Ahorn.History.SelectionSnapshot("Selections", relevantRoom, selections)
    ])

    needsRedraw |= Ahorn.callbackFirstActive(hotkeys, event)

    if haskey(Ahorn.moveDirections, event.keyval)
        needsRedraw |= handleMovement(event)
    end

    if haskey(resizeModifiers, event.keyval)
        needsRedraw |= handleResize(event)
    end

    if haskey(scaleMultipliers, event.keyval) && !Ahorn.modifierControl()
        needsRedraw |= handleScaling(event)
    end

    if event.keyval in addNodeKeys && !Ahorn.modifierControl()
        needsRedraw |= handleAddNodes(event)
    end

    if event.keyval == Ahorn.Gtk.GdkKeySyms.Delete && !Ahorn.modifierControl()
        needsRedraw |= handleDeletion(selections)
    end

    if needsRedraw
        clearDragging!()
        clearResize!()
        mouseMotionAbs(lastX, lastY)

        toolsLayer.redraw = true
        redrawTargetLayer!(targetLayer, layersSelected)

        Ahorn.History.addSnapshot!(snapshot)
    end

    return true
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
# Using names.txt (right click and 'Save Link/Target As...'), a 46K text file
# containing over five-thousand first names, begin by sorting it into
# alphabetical order. Then working out the alphabetical value for each name,
# multiply this value by its alphabetical position in the list to obtain a name
# score.
#
# For example, when the list is sorted into alphabetical order, COLIN, which is
# worth 3 + 15 + 12 + 9 + 14 = 53, is the 938th name in the list. So, COLIN
# would obtain a score of 938 × 53 = 49714.
#
# What is the total of all the name scores in the file?

using ProjectEulerSolutions

# Use sort then map and collect to order names then convert to char array then
# difference from the char 'A'.
function p022_solution(path::String="")::Integer
    all_names = Vector{String}()
    if isfile(path)
        open(path) do file
            for row in eachline(file)
                rownames = split(row, ",")
                append!(all_names, [n[2:end-1] for n in rownames]) # Remove quotes
            end
        end
    end
    sort!(all_names)

    return sum(map(x->sum(collect(x[2]).-'A'.+1)*x[1], enumerate(all_names)))
end

p022 = Problems.Problem(p022_solution)

input_path = "../data/p022_names.txt"

Problems.benchmark(p022, input_path)

ext: jl

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
#################### Monte Carlo Standard Errors ####################

function mcse(x::Vector{T}, method::Symbol=:imse; args...) where {T<:Real}
  method == :bm ? mcse_bm(x; args...) :
  method == :imse ? mcse_imse(x) :
  method == :ipse ? mcse_ipse(x) :
    throw(ArgumentError("unsupported mcse method $method"))
end

function mcse_bm(x::Vector{T}; size::Integer=100) where {T<:Real}
    n = length(x)
    m = div(n, size)
    if m < 2
        throw(
         ArgumentError("iterations are < $(2 * size) and batch size is > $(div(n, 2))")
        )
    end
    mbar = [mean(x[i * size .+ (1:size)]) for i in 0:(m - 1)]
    return sem(mbar)
end

function mcse_imse(x::Vector{T}) where {T<:Real}
    n = length(x)
    m = div(n - 2, 2)
    x_ = map(Float64, x)
    ghat = autocov(x_, [0, 1])
    Ghat = sum(ghat)
    value = -ghat[1] + 2 * Ghat
    for i in 1:m
        Ghat = min(Ghat, sum(autocov(x_, [2 * i, 2 * i + 1])))
        Ghat > 0 || break
        value += 2 * Ghat
    end
    return sqrt(value / n)
end

function mcse_ipse(x::Vector{T}) where {T<:Real}
    n = length(x)
    m = div(n - 2, 2)
    x_ = map(Float64, x)
    ghat = autocov(x_, [0, 1])
    value = ghat[1] + 2 * ghat[2]
    for i in 1:m
        Ghat = sum(autocov(x_, [2 * i, 2 * i + 1]))
        Ghat > 0 || break
        value += 2 * Ghat
    end
    return sqrt(value / n)
end


text:
@enum RuleType nop acc jmp

struct Rule
    instruction::RuleType
    qty::Integer
end

Rule(instruction::RuleType) = Rule(instruction::RuleType, 0)

mutable struct Program
    instructions::Vector{Rule}
    state::Integer
    position::Integer
    iteration_counter::Integer
end

function step_program!(program)
    # we execute the rule we are at now
    current_rule = program.instructions[program.position]

    # here are all the rules - maybe should refactor them and dispatch
    if current_rule.instruction == nop
        program.position += 1
        program.iteration_counter += 1

    elseif current_rule.instruction == acc
        program.state += current_rule.qty
        program.position += 1
        program.iteration_counter += 1

    elseif current_rule.instruction == jmp
        program.position += current_rule.qty
        program.iteration_counter += 1

    else
        throw("Invalid rule $current_rule")
    end
end

function edit_program_until_it_terminates(program::Program)
    for i in 1:length(program.instructions)
        current_rule = program.instructions[i]
        if current_rule.instruction == acc
            continue
        end

        new_program = deepcopy(program)
        if current_rule.instruction == jmp
            new_rule = Rule(nop, current_rule.qty)
            new_program.instructions[i] = new_rule
        elseif current_rule.instruction == nop
            new_rule = Rule(jmp, current_rule.qty)
            new_program.instructions[i] = new_rule
        end

        if run_program_until_duplicate_line(new_program)
            println("Mutation at line $i gets the program to terminate with $(new_program.state)")
        end
    end
end

function run_program_until_duplicate_line(program::Program)
    executed_lines = Set([])
    while !(program.position in executed_lines)
        push!(executed_lines, program.position)
        step_program!(program)
        if program.position == (length(program.instructions) + 1)
            # println("success")
            return true
        end
    end

    # println("Entered loop.  program state: $(program.state) program_iteration $(program.iteration_counter)")
    return false
end

function parse_instruction_into_rule(line)
    # lines look like: nop +0
    instruction_name, qty = split(line, " ")

    qty = parse(Int64, qty)
    if instruction_name == "nop"
        instruction_type = nop
    elseif instruction_name == "acc"
        instruction_type = acc
    elseif instruction_name == "jmp"
        instruction_type = jmp
    else
        throw("Invalid instruction $instruction_name")
    end
    Rule(instruction_type, qty)
end

function build_program(input_file)
    program = open(input_file) do file
        file_string = read(file, String)
        instructions = split(file_string, "\r\n")
        rules = [parse_instruction_into_rule(line) for line in instructions]
        program = Program(rules, 0, 1, 1)
    end
    program
end

p = build_program("day8_input.txt")
edit_program_until_it_terminates(p)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""DEPRECATED 1.9 - This group version of NetworkPolicyPort is deprecated by networking/v1/NetworkPolicyPort.

    IoK8sApiExtensionsV1beta1NetworkPolicyPort(;
        port=nothing,
        protocol=nothing,
    )

    - port::IoK8sApimachineryPkgUtilIntstrIntOrString : If specified, the port on the given protocol.  This can either be a numerical or named port on a pod.  If this field is not provided, this matches all port names and numbers. If present, only traffic on the specified protocol AND port will be matched.
    - protocol::String : Optional.  The protocol (TCP, UDP, or SCTP) which traffic must match. If not specified, this field defaults to TCP.
"""
mutable struct IoK8sApiExtensionsV1beta1NetworkPolicyPort <: SwaggerModel
    port::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: port
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol

    function IoK8sApiExtensionsV1beta1NetworkPolicyPort(;port=nothing, protocol=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyPort, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        validate_property(IoK8sApiExtensionsV1beta1NetworkPolicyPort, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        o
    end
end # type IoK8sApiExtensionsV1beta1NetworkPolicyPort

const _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPort = Dict{Symbol,Symbol}(Symbol("port")=>Symbol("port"), Symbol("protocol")=>Symbol("protocol"))
const _property_types_IoK8sApiExtensionsV1beta1NetworkPolicyPort = Dict{Symbol,String}(Symbol("port")=>"IoK8sApimachineryPkgUtilIntstrIntOrString", Symbol("protocol")=>"String")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPort }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPort))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPort }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiExtensionsV1beta1NetworkPolicyPort[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPort }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1NetworkPolicyPort[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1NetworkPolicyPort)
    true
end

function validate_property(::Type{ IoK8sApiExtensionsV1beta1NetworkPolicyPort }, name::Symbol, val)
end


text:
export CIFAR100
module CIFAR100

using BinDeps

const defdir = joinpath(Pkg.dir("MLDatasets"), "datasets/cifar100")

function getdata(dir)
    mkpath(dir)
    path = download("https://www.cs.toronto.edu/~kriz/cifar-100-binary.tar.gz")
    run(unpack_cmd(path,dir,".gz",".tar"))
end

function readdata(data::Vector{UInt8})
    n = Int(length(data)/3074)
    x = Array(Float64, 3072, n)
    y = Array(Int, 2, n)
    for i = 1:n
        k = (i-1) * 3074 + 1
        y[:,i] = data[k:k+1]
        x[:,i] = data[k+2:k+3073] / 255
    end
    x = reshape(x, 32, 32, 3, n)
    x, y
end

function traindata(dir=defdir)
    file = joinpath(dir, "cifar-100-binary/train.bin")
    isfile(file) || getdata(dir)
    readdata(open(read,file))
end

function testdata(dir=defdir)
    file = joinpath(dir, "cifar-100-binary/test.bin")
    isfile(file) || getdata(dir)
    readdata(open(read,file))
end

end


text:
using DataStructures
using CountMin

N = 100000

function dotest(cms, n::Integer)
    while n > 0
        push!(cms, n-=1)
    end
end

println("Doing CMS")
cms = CountMinSketch{UInt8}(4, 2^28)
dotest(cms, 1000)
@time dotest(cms, N)

println("Doing Counter")
ctr = counter(Int)
dotest(ctr, 1000)
@time dotest(ctr, N)

print("Checking equality... ")
for i in 1:N
    assert(cms[i] == ctr[i])
end
println("All Good!")

function sumtest(cms, n)
    sm = 0
    for i in 1:N
        sm += cms[i]
    end
end

@time sumtest(cms, N)

println("Profiling add/push")
Profile.clear()
@profile dotest(cms, N)
Profile.print(maxdepth=15)

println("Profiling getitem")
Profile.clear()
@profile sumtest(cms, N)
Profile.print(maxdepth=15)


text:
@testset "basic" begin

# test of basic funtionality

#tests for config set/get
cnf = Spark.SparkConf()
cnf["abc"] = "def"
@test cnf["abc"] == "def"

#tests for SparkContext
sc = SparkContext(master="local")
txt = parallelize(sc, ["hello", "world"])
rdd = map_partitions(txt, it -> map(s -> length(s), it))

@test count(rdd) == 2
@test reduce(rdd, +) == 10

close(sc)

end


text:
module PseudoSteppingSchemes

using Statistics
using LinearAlgebra

export ConstantConvergence

using ..EnsembleKalmanInversions: step_parameters
import ..EnsembleKalmanInversions: adaptive_step_parameters

struct ConstantConvergence{T}
    convergence_ratio :: T
end

function volume_ratio(Xⁿ⁺¹, Xⁿ)
    Vⁿ⁺¹ = det(cov(Xⁿ⁺¹, dims=2))
    Vⁿ   = det(cov(Xⁿ,   dims=2))
    return Vⁿ⁺¹ / Vⁿ
end

function adaptive_step_parameters(pseudo_stepping::ConstantConvergence, Xⁿ, Gⁿ, y, Γy, process; Δt=1.0)
    convergence_ratio = pseudo_stepping.convergence_ratio

    # Test step forward
    Xⁿ⁺¹ = step_parameters(Xⁿ, Gⁿ, y, Γy, process; Δt)
    r = volume_ratio(Xⁿ⁺¹, Xⁿ)

    # "Accelerated" fixed point iteration to adjust (pseudo) Δt
    p = 1.1
    iter = 1
    while !isapprox(r, convergence_ratio, atol=0.03, rtol=0.1) && iter < 10
        Δt *= (r / convergence_ratio)^p
        Xⁿ⁺¹ = step_parameters(Xⁿ, Gⁿ, y, Γy, process; Δt)
        r = volume_ratio(Xⁿ⁺¹, Xⁿ)
        iter += 1
    end

    @info "Particles stepped adaptively with convergence ratio $r (target $convergence_ratio)"

    return Xⁿ⁺¹, Δt
end

end # module


text:
struct Framewise{AT} <: ADIAlgorithm
    kernel::AT
    limit
    delta_rot
end

"""
    Framewise(alg; limit=Inf, delta_rot=1)
    Framewise(algs::AbstractVector; limit=Inf, delta_rot=1)

Wrap an algorithm such that the underlying data will be processed frame by frame. For each frame a reference library is created from the data. This reference can be filtered by rejecting frames which have not rotated a sufficient parallactic angle. `delta_rot` sets the required arc length for rotation in units of the FWHM. The number of frames retained can be specified with `limit`, e.g. the 4 closest frames in time with the target frame.

Because the framewise application of an algorithm requires additional information, the following keyword arguments must be provided to [`reconstruct`](@ref) or [`subtract`](@ref).
* `fwhm` - the FWHM of the instrument in pixels. Will be set to the width of a [`MultiAnnulusView`](@ref)
* `angles` - the measured parallactic angles for each frame
* `r` - The radius of the arc to calculate the parallactic angle threshold. Will be set automatically if using [`AnnulusView`](@ref) or [`MultiAnnulusView`](@ref).

In addition, `Framewise` versions of algorithms do not implement [`ADI.fit`](@ref) and do not currently support RDI.

## Annular reduction

In the case of reducing a [`MultiAnnulusView`](@ref), a vector of algorithms can be used, each one corresponding to each annulus of the view. In this case, too, `delta_rot` can be given as a vector or as a tuple. If it is given as a tuple, `delta_rot` will increase linearly from the first value to the last value across each annulus.

# Examples

```julia
julia> cube, angles = # load data

julia> alg = PCA(10) # the algorithm to use on each reference

julia> res = Framewise(alg)(cube, angles; r=10, fwhm=5);

julia> mav = MultiAnnulusView(cube, 5; inner=5);

julia> res_ann = Framewise(alg, delta_rot=(0.1, 1))(mav, angles);
```
"""
Framewise(alg; limit=Inf, delta_rot=1) = Framewise(alg, limit, delta_rot)

function reconstruct(alg::Framewise, cube::AbstractArray{T,3}; angles, fwhm, r, kwargs...) where T
    pa_threshold = compute_pa_thresh(angles, r, fwhm, alg.delta_rot)
    data = flatten(cube)
    S = similar(data)
    @views Threads.@threads for i in axes(data, 1)
        inds = find_angles(angles, i, pa_threshold; limit=alg.limit)
        target = data[i:i, :]
        ref = data[inds, :]
        angs = angles[inds]
        des = fit(alg.kernel, target; kwargs..., ref=ref, angles=angs)
        S[i, :] = reconstruct(des)
    end
    return expand(S)
end

function reconstruct(alg::Framewise, cube::AnnulusView; angles, fwhm, r=_radius(cube), kwargs...)
    pa_threshold = compute_pa_thresh(angles, r, fwhm, alg.delta_rot)
    data = cube(true) # as view
    S = similar(data)
    Threads.@threads for i in axes(S, 1)
        inds = find_angles(angles, i, pa_threshold; limit=alg.limit)
        target = data[i:i, :]
        ref = data[inds, :]
        angs = angles[inds]
        des = fit(alg.kernel, target; kwargs..., ref=ref, angles=angs)
        S[i, :] = reconstruct(des)
    end
    return inverse(cube, S)
end

_radius(cube::AnnulusView) = cube.rmin + (cube.rmax - cube.rmin)/2

function reconstruct(alg::Framewise, cube::MultiAnnulusView; angles, fwhm=cube.width, kwargs...)
    local recons
    anns = eachannulus(cube, true) # as views
    N_ann = length(cube.indices)
    delta_rots = _normalize_deltarot(alg.delta_rot, N_ann)
    @withprogress name="annulus" begin
        i_ann = 0
        recons = map(anns, cube.radii, delta_rots) do ann, r, delta_rot
            pa_threshold = compute_pa_thresh(angles, r, fwhm, delta_rot)
            @debug "PA thresh: $pa_threshold Ann center: $r"
            S = similar(ann)
            @views Threads.@threads for j in axes(S, 1)
                inds = find_angles(angles, j, pa_threshold; limit=alg.limit)
                target = ann[j:j, :]
                ref = ann[inds, :]
                angs = angles[inds]
                des = fit(alg.kernel, target; kwargs..., ref=ref, angles=angs)
                S[j, :] = reconstruct(des)
            end
            # update progress
            i_ann += 1
            @logprogress i_ann/N_ann thresh=pa_threshold
            return S
        end
    end

    return inverse(cube, recons)
end

_normalize_deltarot(delta_rot, N) = Fill(delta_rot, N)
_normalize_deltarot(delta_rot::Tuple, N) = range(delta_rot..., length=N)
_normalize_deltarot(delta_rot::AbstractVector, N) = delta_rot

function reconstruct(alg::Framewise{<:AbstractVector}, cube::MultiAnnulusView; angles, fwhm=cube.width, kwargs...)
    anns = eachannulus(cube, true) # as views
    N_ann = length(cube.indices)
    delta_rots = _normalize_deltarot(alg.delta_rot, N_ann)
    local recons
    @withprogress name="annulus" begin
        i_ann = 0
        recons = map(anns, cube.radii, alg.kernel, delta_rots) do ann, r, _alg, delta_rot
            pa_threshold = compute_pa_thresh(angles, r, fwhm, delta_rot)
            S = similar(ann)
            @views Threads.@threads for j in axes(S, 1)
                @debug "PA thresh: $pa_threshold Ann center: $r"
                inds = find_angles(angles, j, pa_threshold; limit=alg.limit)
                target = ann[j:j, :]
                ref = ann[inds, :]
                angs = angles[inds]
                des = fit(_alg, target; kwargs..., ref=ref, angles=angs)
                S[j, :] = reconstruct(des)
            end
            # update progress
            i_ann += 1
            @logprogress i_ann/N_ann
            return S
        end
    end
    return inverse(cube, recons)
end

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

function compute_pa_thresh(angles, r, fwhm, delta_rot)
    pa_threshold = 2 * atand(delta_rot * fwhm,  2r)
    mid_range = abs(maximum(angles) - minimum(angles)) / 2
    k = mid_range * 0.9
    if pa_threshold ≥ k
        @info "pa threshold $pa_threshold too large, will be set to $k"
        pa_threshold = k
    end
    return pa_threshold
end

function find_angles(angles, idx, thresh; limit=Inf)
    fidx, lidx = firstindex(angles), lastindex(angles)
    p, n = fidx, idx
    for i in fidx:idx-1
        if abs(angles[idx] - angles[i]) < thresh
            p = i
            break
        else
            p += 1
        end
    end
    for k in idx:lidx
        if abs(angles[k] - angles[idx]) > thresh
            n = k
            break
        else
            n += 1
        end
    end
    if isfinite(limit)
        window = limit ÷ 2
        first_half = max(p - window, fidx):p - 1
        last_half = n:min(n + window, lidx) - 1
    else
        first_half = fidx:p - 1
        last_half = n:lidx
    end
    return vcat(first_half, last_half)
end


text:
import Cairo
using Colors
using Base.Test
using Winston
using Compat

module ImageComparisons
    using Winston, Colors
    include("examples.jl")
    include("issues.jl")
    include("plot.jl")
end

function read_png_data(fn::ASCIIString)
    surface = Cairo.read_from_png(fn)
    w = Cairo.width(surface)
    h = Cairo.height(surface)
    p = ccall((:cairo_image_surface_get_data,Cairo._jl_libcairo),
              Ptr{UInt8}, (Ptr{Void},), surface.ptr)
    a = pointer_to_array(convert(Ptr{UInt32},p), (convert(Int,w),convert(Int,h)))
    copy(a)
end

function img_dist{T<:Array{UInt32,2}}(img1::T, img2::T)
    @assert size(img1) == size(img2)
    s = 0.
    for i = 1:length(img1)
        a = convert(RGB, RGB24(img1[i]))
        b = convert(RGB, RGB24(img2[i]))
        s += colordiff(a, b)
    end
    s
end

function main()
    root = dirname(task_local_storage()[:SOURCE_PATH])
    dir1,dir2 = "$root/_baseline","$root/_current"
    isdir(dir1) || mkdir(dir1)
    isdir(dir2) || mkdir(dir2)

    for name in sort(names(ImageComparisons))

        name == :ImageComparisons && continue
        p = eval(:((ImageComparisons.$name)()))

        fn1 = "$dir1/$name.png"
        fn2 = "$dir2/$name.png"
        isfile(fn1) || savefig(p, fn1)
        savefig(p, fn2)

        img1 = read_png_data(fn1)
        img2 = read_png_data(fn2)

        d = img_dist(img1, img2)
        println("$name $d")
        d < 0.1 || error("$name failed")
    end
end

main()


text:
# restart processes
if nprocs() > 1
    rmprocs(workers()) # remove all worker processes
end
wpids = addprocs(2) # add processes

println("Spawned ", nprocs(), " processes, ", nworkers()," workers")
println("Proc IDs: ", procs())

# load LMDB on all processes
@everywhere using LMDB

# create a sample database
nsamples = 100
dbname = "simpleseq.db"
isdir(dbname) && rm(dbname, recursive=true)
!isdir(dbname) && mkdir(dbname)

# the data are just {1:1, 2:2 ... }
environment(dbname) do env
    start(env) do txn
        open(txn) do dbi
            for i=1:nsamples
                put!(txn, dbi, string(i), string(i))
            end
        end
        commit(txn)
    end
end

# helper functions
@everywhere function getSamplesFromDb(env, idxs::Array{Int})
    txn = start(env)
    dbi = open(txn)
    xs = Int[]
    println(idxs)
    for idx in idxs
        key = string(idx)
        val = get(txn, dbi, key, String);
        println("k:$key, v:$val")
        val = parse(Int, val)
        push!(xs, val)
    end
    abort(txn)
    close(env, dbi)
    return xs
end

@everywhere function miniBatchSum(idxs, dbname::String)
    # open the database
    println("Opening ", dbname)
    xs = environment(dbname) do env
        getSamplesFromDb(env, idxs)
    end

    # the cost is the sum of all the values we get from the db
    cost = 0.0
    for x in xs
        cost += x;
    end

    return cost
end

# the following single process call works
println("miniBatchSum([1,2,3], dbname)")
println(miniBatchSum([1,2,3], dbname))

# the following (which does it in parallel) does not work
# we generate some ids to split across the nodes
# each node will process sample_size values
# the ids are put into proc_idxs
sample_size = 10;
idxs = randperm(nsamples);
idxs = idxs[1:(nworkers()*sample_size)]
proc_idxs = Any[]
st_idx = 1;
en_idx = sample_size;
for i=1:nworkers()
    push!(proc_idxs, idxs[st_idx:en_idx]);
    st_idx = en_idx+1;
    en_idx = en_idx+sample_size;
end
println(proc_idxs)

# spawn and run across all worker nodes
k = 1;
remrefs = Array{Future}(nworkers());
for proc in workers()
    println("Remote call to: ", proc);
    remrefs[k] = remotecall(miniBatchSum, proc, proc_idxs[k], dbname);
    k += 1;
end

# collect the results
k = 1;
results = Array{Float64}(nworkers());
for k = 1:length(remrefs)
    wait(remrefs[k]);
    results[k] = fetch(remrefs[k]);
    k += 1;
end
println(results)


text:
# run all 2 processor tests in single session

args_orig = copy(ARGS)

cd("./advection")
include(joinpath(pwd(), "runtests_parallel2.jl"))

cd("../euler")
resize!(ARGS, length(args_orig))  # reset to default
copy!(ARGS, args_orig)
include(joinpath(pwd(), "runtests_parallel2.jl"))

cd("../")


text:
# Tests on Wishart distributions

using Distributions
using Base.Test

V = [[2. 1.], [1. 2.]]
W = Wishart(3., V)

# logdet

@test_approx_eq expected_logdet(W) 1.9441809588650447

# entropy

@test_approx_eq entropy(W) 7.178942679971454


text:
@reexport module CaretLearners

export CaretLearner,fit!,transform!

using TSML

import TSML.TSMLTypes.fit!   # importing to overload
import TSML.TSMLTypes.transform! # importing to overload

using RCall

function initlibs()
    #packages = ["caret","e1071","gam","randomForest","nnet","kernlab","grid","MASS","pls"]
    packages = ["caret","randomForest"]
    for pk in packages
        rcall(:library,pk,"lib=.libPaths()")
    end
end

mutable struct CaretLearner <: TSLearner
    model
    args

    function CaretLearner(args=Dict())
        #fitControl=:(R"trainControl(method = 'repeatedcv',number = 5)")
        fitControl="trainControl(method = 'none')"
        default_args = Dict(
            :output => :class,
            :learner => "rf",
            :fitControl=>fitControl,
            :impl_args => Dict()
        )
        initlibs()
        new(nothing,mergedict(default_args,args))
    end
end

function fit!(crt::CaretLearner,x::T,y::Vector) where  {T<:Union{Vector,Matrix,DataFrame}}
    xx = x |> DataFrame
    yy = y |> Vector
    rres = rcall(:train,xx,yy,method=crt.args[:learner],trControl = reval(crt.args[:fitControl]))
    #crt.model = R"$rres$finalModel"
    crt.model = rres
end

function transform!(crt::CaretLearner,x::T) where  {T<:Union{Vector,Matrix,DataFrame}}
    xx = x |> DataFrame
    res = rcall(:predict,crt.model,xx) #in robj
    return rcopy(res) # return extracted robj
end

function caretrun()
    crt = CaretLearner(Dict(:learner=>"rf") )
    iris=getiris()
    x=iris[:,1:4]  |> Matrix
    y=iris[:,5] |> Vector
    fit!(crt,x,y)
    print(crt.model)
    transform!(crt,x) |> collect
end

end


text:
using ZigZagBoomerang
using Test
using ZigZagBoomerang: poisson_time
"""
int (max(a + b*t, 0) + c) dt from 0 to pi
"""
function F((a,b,c), π)
    if a <= 0 && (a + π*b <= 0)
        return π*c
    elseif a > 0 && (a + π*b < 0)
        return π*c - a^2/(2*b)
    elseif a > 0 && (a + π*b >= 0)
        return 1/2*π*(2*a + π*b + 2*c)
    else # a <= 0 && (a + π*b > 0)
        return a^2/(2*b) + π*a + (π^2*b)/2 + π*c
    end
end
@testset "poisson" begin
    for k in 1:100
        a, b = 2rand(2) .- 1
        u = rand()
        s = poisson_time((a,b), u)
        if s == Inf
            t = F((a,b,false), s)
            @test t < -log(u)
        else
            t = F((a,b,false), s)
            @test t ≈ -log(u)
        end
    end
end

@testset "poisson" begin
    for k in 1:100
        a, b = 2rand(2) .- 1
        c = rand()
        u = rand()
        s = poisson_time((a,b,c), u)
        if s == Inf
            t = F((a,b,c), s)
            @test t < -log(u)
        else

            t = F((a,b,c), s)
            @test t ≈ -log(u)
        end
    end
end


text:
# Unit tests for api/md_simulator.jl

struct MockSystem <: AbstractSystem{3} end
struct MockSimulator <: MolecularDynamicsSimulator end
struct MockPotential <: InteratomicPotentials.ArbitraryPotential end

@testset "api/md_simulator.jl" begin
    @test_throws Atomistic.UnimplementedError{MockSimulator} simulate(MockSystem(), MockSimulator(), MockPotential())
end


text:
for N in [Float64, Float32, Rational{Int}]
    c1 = N[1, 2]
    c2 = N[3, 4, 5]
    Z1 = Zonotope(c1, N[1 0; 0 1])
    Z2 = Zonotope(c2, N[1 0 0; 0 1 0; 0 0 1])
    H1 = Hyperrectangle(c1, N[1, 1])
    H2 = Hyperrectangle(c2, N[1, 1, 1])
    B1 = BallInf(c1, N(1))
    B2 = BallInf(c2, N(1))

    Z1Z2 = cartesian_product(Z1, Z2)
    Z1H2 = cartesian_product(Z1, H2)
    H1Z2 = cartesian_product(H1, Z2)
    H1H2 = cartesian_product(H1, H2)
    H1B2 = cartesian_product(H1, B2)
    B1B2 = cartesian_product(B1, B2)
    for X in [Z1Z2, Z1H2, H1Z2]
        @test X isa Zonotope{N}
    end
    for X in [H1H2, H1B2, B1B2]
        @test X isa Hyperrectangle{N}
    end
    if test_suite_polyhedra || N <: AbstractFloat
        for X in [Z1Z2, Z1H2, H1Z2, H1H2, H1B2]
            @test isequivalent(X, B1B2)
        end
    end

    # Universe is handled correctly (#2735)
    P = HalfSpace(N[1, 2], N(3))
    U = Universe{N}(2)
    @test isequivalent(cartesian_product(P, U), HalfSpace(N[1, 2, 0, 0], N(3)))
    @test isequivalent(cartesian_product(U, P), HalfSpace(N[0, 0, 1, 2], N(3)))
    @test cartesian_product(U, U) == Universe{N}(4)

    # hyperrectangle and Universe
    H = Hyperrectangle(N[1, 2], N[3, 4])
    HU = cartesian_product(H, U)
    @test HU isa HPolyhedron{N, SingleEntryVector{N}} && isequivalent(HU,
        HPolyhedron([HalfSpace(N[1, 0, 0, 0], N(4)),
                     HalfSpace(N[-1, 0, 0, 0], N(2)),
                     HalfSpace(N[0, 1, 0, 0], N(6)),
                     HalfSpace(N[0, -1, 0, 0], N(2))]))
    UH = cartesian_product(U, H)
    @test UH isa HPolyhedron{N, SingleEntryVector{N}} && isequivalent(UH,
        HPolyhedron([HalfSpace(N[0, 0, 1, 0], N(4)),
                     HalfSpace(N[0, 0, -1, 0], N(2)),
                     HalfSpace(N[0, 0, 0, 1], N(6)),
                     HalfSpace(N[0, 0, 0, -1], N(2))]))

    # HalfSpace and Universe
    H1 = HalfSpace(N[0, 4, 0], N(5))
    H2 = HalfSpace(sparsevec([2], N[4], 3), N(5))
    H3 = HalfSpace(SingleEntryVector(2, 3, N(4)), N(5))
    @test isequivalent(H1, H2) && isequivalent(H1, H3)
    H1a = cartesian_product(H1, Universe{N}(2))
    H1p = cartesian_product(Universe{N}(2), H1)
    @test isequivalent(H1a, cartesian_product(H2, Universe{N}(2))) &&
          isequivalent(H1a, cartesian_product(H3, Universe{N}(2))) &&
          isequivalent(H1a, HalfSpace(N[0, 4, 0, 0, 0], N(5)))
    @test isequivalent(H1p, cartesian_product(Universe{N}(2), H2)) &&
          isequivalent(H1p, cartesian_product(Universe{N}(2), H3)) &&
          isequivalent(H1p, HalfSpace(N[0, 0, 0, 4, 0], N(5)))

    # polyhedron and Universe
    P = HPolyhedron([HalfSpace(N[1], N(2))])
    PU = cartesian_product(P, U)
    @test isequivalent(PU, HPolyhedron([HalfSpace(N[1, 0, 0], N(2))]))
    UP = cartesian_product(U, P)
    @test isequivalent(UP, HPolyhedron([HalfSpace(N[0, 0, 1], N(2))]))
end


text:
## define variables...
INSTALL_GPU = Sys.islinux() && haskey(ENV, "GPU") && ENV["GPU"] in [1, "1"]

CONDA = ""
if Sys.iswindows()
    CONDA = "$(JULIA_ADCME_DIR)/.julia/adcme/Scripts/conda.exe"
else
    CONDA = "$(JULIA_ADCME_DIR)/.julia/adcme/bin/conda"
end

INSTALLER = ""
if Sys.islinux()
    INSTALLER = "Miniconda3-py37_4.8.3-Linux-x86_64.sh"
elseif Sys.isapple()
    INSTALLER = "Miniconda3-py37_4.8.3-MacOSX-x86_64.sh"
else
    INSTALLER = "Miniconda3-py37_4.8.3-Windows-x86_64.exe"
end

function check_install()
    if Sys.iswindows()
    else
        for bin in ["unzip", "ninja"]
            if !isfile(joinpath("$(JULIA_ADCME_DIR)/.julia/adcme/bin"), bin)
                return false
            end
        end
    end
    return true
end

function install_conda()
    PWD = pwd()
    cd("$(JULIA_ADCME_DIR)/.julia/")
    if !(INSTALLER in readdir("."))
        @info "Downloading miniconda installer..."
        download("https://repo.anaconda.com/miniconda/"*INSTALLER, INSTALLER)
    end

    if isdir("adcme")
        if FORCE_REINSTALL_ADCME
            error("""ADCME dependencies directory already exist, and you indicate FORCE_REINSTALL_ADCME=true.
    Please (1) quit Julia, (2) remove the path $(joinpath(pwd(), "adcme")) and (3) rebuild ADCME.""")
        else
            @info "ADCME dependencies have already been installed."
            return
        end
    end

    @info "Installing miniconda..."
    if Sys.iswindows()
        run(`cmd /c start /wait "" $INSTALLER /InstallationType=JustMe /RegisterPython=0 /S /D=$(JULIA_ADCME_DIR)\\.julia\\adcme`)
    else
        run(`bash $INSTALLER -f -b -p adcme`)
    end
    cd(PWD)
end

function install_conda_envs()
    cd(@__DIR__)

    ENV_ = copy(ENV)
    if Sys.iswindows()
        platform = "windows"
        ENV_["PATH"] = "$(JULIA_ADCME_DIR)/.julia/adcme/Scripts;$(JULIA_ADCME_DIR)/.julia/adcme/Library/bin;$(JULIA_ADCME_DIR)/.julia/adcme/" * ENV_["PATH"]
    elseif Sys.islinux()
        platform = INSTALL_GPU ? "linux-gpu" : "linux"
    else
        platform = "osx"
    end

    if check_install()
        if (platform in ["windows", "linux", "osx"] && occursin("tensorflow", read(`$CONDA list`, String)))
            return
        end

        if platform == "linux-gpu"
            if occursin("tensorflow-gpu", read(`$CONDA list`, String))
                return
            elseif occursin("tensorflow", read(`$CONDA list`, String))
                @error """You have already installed tensorflow-cpu. In order to install tensorflow-gpu, try the following steps:
- Quit Julia and remove the directory

$(JULIA_ADCME_DIR)/.julia/adcme

- In Julia, run the following command

```
ENV["GPU"] = 1
using Pkg; Pkg.build("ADCME")
```
"""
            end
        end

    end
    @info "Installing conda dependencies..."
    run(setenv(`$CONDA env update -n base --file $platform.yml`, ENV_))
end

install_conda()
install_conda_envs()


text:
function pivot_longer(df::AbstractDataFrame,bycol::Symbol,cols::Vector{Symbol};period = :Period, separator = true)
    # temp = nothing
     df_longer = nothing
     m = nothing
    @inbounds for col ∈ cols
        try
        temp = select(df,Regex(String(col)))
        @inbounds for colname ∈ names(temp)
          if separator
              m = match(Regex("($(col))(?<Value>\\d+)"),String(colname))
          else
              m = match(Regex("($(col))(?<Value>\\d+)"),String(colname))
          end
          if m == nothing
            println("Column $(col) has no regex match")
          else
            rename!(temp,colname=>Symbol(m[2]))
          end
        end
        temp[!,bycol] = df[!,bycol]
        temp = stack(temp,Not(bycol),variable_eltype=String)
        rename!(temp,:variable => period, :value => Symbol(m[1]))
        if df_longer == nothing
            df_longer = temp
        else
            df_longer = outerjoin(df_longer, temp, on = [bycol,period],makeunique=true)
        end
        catch
            @warn "Problem with column $col"
        end
    end
    return df_longer
end

match(r"(age)(?<Value>\d+)",String(:age67))

function reshape_results!(mc,D)
  N = D.N
  J = D.J
  temp = zeros(D.N*D.J)
  Data = DataFrame(temp = temp)
  for i in fieldnames(typeof(mc))
   Data[!,i] = vec(getfield(mc,i))
  end
  select!(Data,Not(:temp))
  return Data
end

function print_struct(S)
  s = maximum(length.(string.(fieldnames(typeof(S)))))
  for i in fieldnames(typeof(S))
  @info "$(rpad(i,s," ")) = $(round.(getfield(S,i),digits=2))"
  end
end

function print_parameters(P)
  for i ∈ fieldnames(typeof(P))
    println("--------------------")
    @info typeof(getfield(P,i))
    println("--------------------")
    print_struct(getfield(P,i))
  end
end

"""
```
percentiles!(data,by_variable,variable;pctls = [0.05, 0.5, 0.95])
```

Example:

data = percentiles!(df,:age,:income,pctls = [0.1,0.9])

Returns additional columns:

data[:income_p1, :income_p2]

data[:bin_income] # with values 1, 2, 3

"""
function percentiles!(data,by_variable,variable;pctls = [0.05, 0.5, 0.95])
  if length(pctls) == 2
    temp = combine(groupby(data,by_variable), variable => (x -> [quantile(x, pctls)]) => [:pctl1, :pctl2])
    data = innerjoin(data,temp,on = by_variable)
    name = Symbol("bin_$variable")
    data[!,name].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] <= data[i,:pctl1]
         data[i,name] = 1
      elseif data[i,variable]> data[i,:pctl1] && data[i,variable]<= data[i,:pctl2]
         data[i,name] = 2
      elseif data[i,variable] > data[i,:pctl2]
         data[i,name] = 3
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_p1"),:pctl2 => Symbol("$(variable)_p2"))
  elseif length(pctls) == 3
     temp = combine(groupby(data,by_variable), variable => (x -> [quantile(x, pctls)]) => [:pctl1, :pctl2, :pctl3])
     data = innerjoin(data,temp,on = by_variable)
     name = Symbol("bin_$variable")
     data[!,name].= 0
     for i ∈ eachindex(data[!,variable])
       if data[i,variable] <= data[i,:pctl1]
          data[i,name] = 1
       elseif data[i,variable]> data[i,:pctl1] && data[i,variable]<= data[i,:pctl2]
          data[i,name] = 2
       elseif data[i,variable]> data[i,:pctl2] && data[i,variable]<= data[i,:pctl3]
          data[i,name] = 3
       elseif data[i,variable] > data[i,:pctl3]
          data[i,name] = 4
       end
     end
     rename!(data,:pctl1 => Symbol("$(variable)_p1"),:pctl2 => Symbol("$(variable)_p2"),:pctl3 => Symbol("$(variable)_p3"))
  elseif length(pctls) == 4
    temp = combine(groupby(data,by_variable), variable => (x -> [quantile(x, pctls)]) => [:pctl1, :pctl2, :pctl3, :pctl4])
    data = innerjoin(data,temp,on = by_variable)
    name = Symbol("bin_$variable")
    data[!,name].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] <= data[i,:pctl1]
         data[i,name] = 1
      elseif data[i,variable]> data[i,:pctl1] && data[i,variable]<= data[i,:pctl2]
         data[i,name] = 2
      elseif data[i,variable]> data[i,:pctl2] && data[i,variable]<= data[i,:pctl3]
         data[i,name] = 3
      elseif data[i,variable] > data[i,:pctl3]  && data[i,variable]<= data[i,:pctl4]
         data[i,name] = 4
      elseif data[i,variable] > data[i,:pctl4]
        data[i,name] = 5
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_p1"),:pctl2 => Symbol("$(variable)_p2"),:pctl3 => Symbol("$(variable)_p3"),:pctl4 => Symbol("$(variable)_p4"))
  elseif length(pctls) == 9
    temp = combine(groupby(data,by_variable), variable => (x -> [quantile(x, pctls)]) => [:pctl1, :pctl2, :pctl3, :pctl4, :pctl5, :pctl6, :pctl7, :pctl8, :pctl9])

    data = innerjoin(data,temp,on = by_variable)
    name = Symbol("bin_$variable")
    data[!,name].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] <= data[i,:pctl1]
        data[i,name] = 1
      elseif data[i,variable]> data[i,:pctl1] && data[i,variable]<= data[i,:pctl2]
        data[i,name] = 2
      elseif data[i,variable]> data[i,:pctl2] && data[i,variable]<= data[i,:pctl3]
        data[i,name] = 3
      elseif data[i,variable] > data[i,:pctl3]  && data[i,variable]<= data[i,:pctl4]
        data[i,name] = 4
      elseif data[i,variable] > data[i,:pctl4]  && data[i,variable]<= data[i,:pctl5]
        data[i,name] = 5
      elseif data[i,variable] > data[i,:pctl5]  && data[i,variable]<= data[i,:pctl6]
        data[i,name] = 6
      elseif data[i,variable] > data[i,:pctl6]  && data[i,variable]<= data[i,:pctl7]
        data[i,name] = 7
      elseif data[i,variable] > data[i,:pctl7]  && data[i,variable]<= data[i,:pctl8]
        data[i,name] = 8
      elseif data[i,variable] > data[i,:pctl8]  && data[i,variable]<= data[i,:pctl9]
        data[i,name] = 9
      elseif data[i,variable] > data[i,:pctl9]
        data[i,name] = 10
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_p1"),:pctl2 => Symbol("$(variable)_p2"),:pctl3 => Symbol("$(variable)_p3"),:pctl4 => Symbol("$(variable)_p4"),:pctl5 => Symbol("$(variable)_p5"),
    :pctl6 => Symbol("$(variable)_p6"),:pctl7 => Symbol("$(variable)_p7"),:pctl8 => Symbol("$(variable)_p8"),:pctl9 => Symbol("$(variable)_p9"))
    else
    error("Function supports 2, 3, 4 and 9 number of percentiles")
  end
  return data
end

function percentiles!(data,variable;pctls = [0.05, 0.5, 0.95])
  if length(pctls) == 2
      for i ∈ 1: length(pctls)
          data[!,Symbol("pctl$i")] .= quantile(data[!,variable],pctls[i])
      end
      name = Symbol("bin_$variable")
      data[!,Symbol("bin_$variable")].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] <= data[i,:pctl1]
         data[i,name] = 1
      elseif data[i,variable]> data[i,:pctl1] && data[i,variable]<= data[i,:pctl2]
         data[i,name] = 2
      elseif data[i,variable] > data[i,:pctl2]
         data[i,name] = 3
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_p1"),:pctl2 => Symbol("$(variable)_p2"))
  elseif length(pctls) == 3
      for i ∈ 1: length(pctls)
          data[!,Symbol("pctl$i")] .= quantile(data[!,variable],pctls[i])
      end
      name = Symbol("bin_$variable")
      data[!,Symbol("bin_$variable")].= 0
     for i ∈ eachindex(data[!,variable])
       if data[i,variable] <= data[i,:pctl1]
          data[i,name] = 1
       elseif data[i,variable]> data[i,:pctl1] && data[i,variable]<= data[i,:pctl2]
          data[i,name] = 2
       elseif data[i,variable]> data[i,:pctl2] && data[i,variable]<= data[i,:pctl3]
          data[i,name] = 3
       elseif data[i,variable] > data[i,:pctl3]
          data[i,name] = 4
       end
     end
     rename!(data,:pctl1 => Symbol("$(variable)_p1"),:pctl2 => Symbol("$(variable)_p2"),:pctl3 => Symbol("$(variable)_p3"))
  elseif length(pctls) == 4
      for i ∈ 1: length(pctls)
          data[!,Symbol("pctl$i")] .= quantile(data[!,variable],pctls[i])
      end
      name = Symbol("bin_$variable")
      data[!,Symbol("bin_$variable")].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] <= data[i,:pctl1]
         data[i,name] = 1
      elseif data[i,variable]> data[i,:pctl1] && data[i,variable]<= data[i,:pctl2]
         data[i,name] = 2
      elseif data[i,variable]> data[i,:pctl2] && data[i,variable]<= data[i,:pctl3]
         data[i,name] = 3
      elseif data[i,variable] > data[i,:pctl3]  && data[i,variable]<= data[i,:pctl4]
         data[i,name] = 4
      elseif data[i,variable] > data[i,:pctl4]
        data[i,name] = 5
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_p1"),:pctl2 => Symbol("$(variable)_p2"),:pctl3 => Symbol("$(variable)_p3"),:pctl4 => Symbol("$(variable)_p4"))
  elseif length(pctls) == 9
      for i ∈ 1: length(pctls)
          data[!,Symbol("pctl$i")] .= quantile(data[!,variable],pctls[i])
      end
      name = Symbol("bin_$variable")
      data[!,Symbol("bin_$variable")].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] <= data[i,:pctl1]
        data[i,name] = 1
      elseif data[i,variable]> data[i,:pctl1] && data[i,variable]<= data[i,:pctl2]
        data[i,name] = 2
      elseif data[i,variable]> data[i,:pctl2] && data[i,variable]<= data[i,:pctl3]
        data[i,name] = 3
      elseif data[i,variable] > data[i,:pctl3]  && data[i,variable]<= data[i,:pctl4]
        data[i,name] = 4
      elseif data[i,variable] > data[i,:pctl4]  && data[i,variable]<= data[i,:pctl5]
        data[i,name] = 5
      elseif data[i,variable] > data[i,:pctl5]  && data[i,variable]<= data[i,:pctl6]
        data[i,name] = 6
      elseif data[i,variable] > data[i,:pctl6]  && data[i,variable]<= data[i,:pctl7]
        data[i,name] = 7
      elseif data[i,variable] > data[i,:pctl7]  && data[i,variable]<= data[i,:pctl8]
        data[i,name] = 8
      elseif data[i,variable] > data[i,:pctl8]  && data[i,variable]<= data[i,:pctl9]
        data[i,name] = 9
      elseif data[i,variable] > data[i,:pctl9]
        data[i,name] = 10
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_p1"),:pctl2 => Symbol("$(variable)_p2"),:pctl3 => Symbol("$(variable)_p3"),:pctl4 => Symbol("$(variable)_p4"),:pctl5 => Symbol("$(variable)_p5"),
    :pctl6 => Symbol("$(variable)_p6"),:pctl7 => Symbol("$(variable)_p7"),:pctl8 => Symbol("$(variable)_p8"),:pctl9 => Symbol("$(variable)_p9"))
    else
    error("Function supports 2, 3, 4 and 9 number of percentiles")
  end
  return data
end

function at_percentiles!(data,by_variable,variable;pctls = [0.05, 0.5, 0.95])
  if length(pctls) == 2
    temp = by(data, by_variable) do df
       (pctl1 = quantile(df[!,variable],pctls[1]),pctl2 = quantile(df[!,variable],pctls[2]))
    end
    data = join(data,temp,on = by_variable)
    name = Symbol("bin_at_$variable")
    data[!,name].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] ≈ data[i,:pctl1]
         data[i,name] = 1
      elseif data[i,variable]≈ data[i,:pctl2]
         data[i,name] = 2
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_at_p1"),:pctl2 => Symbol("$(variable)_at_p2"))
  elseif length(pctls) == 3
     temp = by(data, by_variable) do df
        (pctl1 = quantile(df[!,variable],pctls[1]),pctl2 = quantile(df[!,variable],pctls[2]), pctl3 = quantile(df[!,variable],pctls[3]))
     end
     data = join(data,temp,on = by_variable)
     name = Symbol("bin_at_$variable")
     data[!,name].= 0
     for i ∈ eachindex(data[!,variable])
       if data[i,variable] ≈ data[i,:pctl1]
          data[i,name] = 1
       elseif  data[i,variable] ≈ data[i,:pctl2]
          data[i,name] = 2
       elseif data[i,variable]≈ data[i,:pctl3]
          data[i,name] = 3
       end
     end
     rename!(data,:pctl1 => Symbol("$(variable)_at_p1"),:pctl2 => Symbol("$(variable)_at_p2"),:pctl3 => Symbol("$(variable)_at_p3"))
  elseif length(pctls) == 4
    temp = by(data, by_variable) do df
       (pctl1 = quantile(df[!,variable],pctls[1]),pctl2 = quantile(df[!,variable],pctls[2]), pctl3 = quantile(df[!,variable],pctls[3]), pctl4 = quantile(df[!,variable],pctls[4]))
    end
    data = join(data,temp,on = by_variable)
    name = Symbol("bin_at_$variable")
    data[!,name].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] ≈ data[i,:pctl1]
         data[i,name] = 1
      elseif data[i,variable] ≈ data[i,:pctl2]
         data[i,name] = 2
      elseif data[i,variable]≈ data[i,:pctl3]
         data[i,name] = 3
      elseif data[i,variable] ≈ data[i,:pctl4]
         data[i,name] = 4
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_at_p1"),:pctl2 => Symbol("$(variable)_at_p2"),:pctl3 => Symbol("$(variable)_at_p3"),:pctl4 => Symbol("$(variable)_at_p4"))
  elseif length(pctls) == 9
    temp = by(data, by_variable) do df
       (pctl1 = quantile(df[!,variable],pctls[1]),pctl2 = quantile(df[!,variable],pctls[2]), pctl3 = quantile(df[!,variable],pctls[3]), pctl4 = quantile(df[!,variable],pctls[4]), pctl5 = quantile(df[!,variable],pctls[5]),
       pctl6 = quantile(df[!,variable],pctls[6]),pctl7 = quantile(df[!,variable],pctls[7]), pctl8 = quantile(df[!,variable],pctls[8]), pctl9 = quantile(df[!,variable],pctls[9]))
    end
    data = join(data,temp,on = by_variable)
    name = Symbol("bin_at_$variable")
    data[!,name].= 0
    for i ∈ eachindex(data[!,variable])
      if data[i,variable] ≈ data[i,:pctl1]
        data[i,name] = 1
      elseif data[i,variable] ≈ data[i,:pctl2]
        data[i,name] = 2
      elseif data[i,variable] ≈ data[i,:pctl3]
        data[i,name] = 3
      elseif data[i,variable] ≈ data[i,:pctl4]
        data[i,name] = 4
      elseif data[i,variable] ≈ data[i,:pctl5]
        data[i,name] = 5
      elseif data[i,variable] ≈ data[i,:pctl6]
        data[i,name] = 6
      elseif data[i,variable] ≈ data[i,:pctl7]
        data[i,name] = 7
      elseif data[i,variable] ≈ data[i,:pctl8]
        data[i,name] = 8
      elseif data[i,variable] ≈ data[i,:pctl9]
        data[i,name] = 9
      end
    end
    rename!(data,:pctl1 => Symbol("$(variable)_at_p1"),:pctl2 => Symbol("$(variable)_at_p2"),:pctl3 => Symbol("$(variable)_at_p3"),:pctl4 => Symbol("$(variable)_at_p4"),:pctl5 => Symbol("$(variable)_at_p5"),
    :pctl6 => Symbol("$(variable)_at_p6"),:pctl7 => Symbol("$(variable)_at_p7"),:pctl8 => Symbol("$(variable)_at_p8"),:pctl9 => Symbol("$(variable)_at_p9"))
    else
    error("Function supports 2, 3, 4 and 9 number of percentiles")
  end
  return data
end

function tab(x)
  x = collect(skipmissing(x))
  temp = DataFrame(Variable = collect(extrema(unique(x))[1]:extrema(unique(x))[2]), Proportion = proportions(x,extrema(unique(x))[1]:extrema(unique(x))[2]))
  temp = temp[temp.Proportion.!=0.0,:]
  return temp
end

function tab(x,wt)
    df = DataFrame(x1 = x, wt1 = wt)
    dropmissing!(df)
    df.wt1 = weights(df.wt1)
  temp = DataFrame(Variable = collect(extrema(unique(df.x1))[1]:extrema(unique(df.x1))[2]), Proportion = proportions(df.x1,extrema(unique(df.x1))[1]:extrema(unique(df.x1))[2],df.wt1))
  temp = temp[temp.Proportion.!=0.0,:]
  return temp
end

function writeGzip(df,path)
  open(GzipCompressorStream, path, "w") do stream
      CSV.write(stream, df)
  end
end

function readGzip(path;head = 1)
  df = CSV.File(transcode(GzipDecompressor, Mmap.mmap(path)),header = head) |> DataFrame
end

function createAgeGroups(df::DataFrame,agemin::Int64,agemax::Int64,step::Int64,agevar::Symbol,varname::Symbol)
    agerange = collect(agemin:step:agemax)
    df[!,varname] .= missing
    N = length(agerange)
    for i in 1:nrow(df)
        for j in 1:N-1
            if !ismissing(df[i,agevar])
                if df[i,agevar] ≥ agerange[j]  && df[i,agevar] < agerange[j+1]
                    df[i,varname] = j
                end
            end
        end
    end
    return df
end


text:
#
# Run this file to generate Figure 6 in
#
#   D. Neumann, T. Wiese, and W. Utschick, Learning the MMSE Channel Estimator,
#   IEEE Transactions on Signal Processing, 2018.
#

push!(LOAD_PATH,".")
using DataFrames
using CSV
import SCM3GPP; const scm = SCM3GPP
import CondNormalMMSE; const mmse = CondNormalMMSE
import CondNormalTF; const cntf = CondNormalTF

include("sim_helpers.jl")
include("more_estimators.jl") # OMP and Genie MMSE

verbose = true
#-------------------------------------
# Simulation parameters
#
write_file = true
filename   = "results/figure6.csv"
nBatches   = 100
nBatchSize = 100

#-------------------------------------
# Channel Model
#
snr        = 0 # [dB]
antennas   = [8,16,32,64,96,128]
AS         = 2.0 # standard deviation of Laplacian (angular spread)
nCoherence = 1
Channel    = scm.SCMMulti(pathAS=AS, nPaths=3)
# method that generates "nBatches" channel realizations
get_channel(nAntennas, nCoherence, nBatches) = scm.generate_channel(Channel, nAntennas, nCoherence=nCoherence, nBatches = nBatches)
# method that samples C_delta from delta prior
get_cov(nAntennas) = scm.toeplitzHe( scm.generate_channel(Channel, nAntennas, nCoherence=1)[2][:] )
# get circulant vector that generates all covariance matrices for arbitrary delta (here: use delta=0)
get_circ_cov_generator(nAntennas) = real(scm.best_circulant_approximation(scm.scm_channel([0.0],[1.0],nAntennas,AS=AS)[2]))

#-------------------------------------
# Learning Algorithm parameters
#
learning_rates_relu    = 1e-4*64./antennas # make learning rates dependend on nAntennas
learning_rates_softmax = 1e-3*ones(antennas)
nLayers = 2
nLearningBatches   = 6000
nLearningBatchSize = 20

results = DataFrame()
# read results from previous run
if isfile(filename)
    results = CSV.read(filename)
end

srand(size(results,1)) # use number of entries in results as seed for random number generator
nn_est = Dict{Symbol,Any}()
for iAntenna in 1:length(antennas)
    nAntennas     = antennas[iAntenna]

    verbose && println("Simulating with ", nAntennas, " antennas")

    # Conditionally normal estimators
    cn_est = Dict{Symbol,Any}()
    cn_est[:FastMMSE] = mmse.FastMMSE(snr, get_circ_cov_generator(nAntennas))
    cn_est[:CircMMSE] = mmse.StructuredMMSE(snr, () -> get_cov(nAntennas), nSamples=16*nAntennas, transform = circ_trans)
    cn_est[:ToepMMSE] = mmse.StructuredMMSE(snr, () -> get_cov(nAntennas), nSamples=16*nAntennas, transform = toep_trans)
    cn_est[:CircML]   = mmse.MLEst(snr, transform = circ_trans)

    # Network estimators
    if iAntenna == 1
        nn_est[:CircReLU] = cntf.ConvNN(nLayers, nAntennas, transform = circ_trans, learning_rate = learning_rates_relu[iAntenna])
        nn_est[:ToepReLU] = cntf.ConvNN(nLayers, nAntennas, transform = toep_trans, learning_rate = learning_rates_relu[iAntenna])
        nn_est[:CircSoftmax] = cntf.ConvNN(nLayers, nAntennas, transform = circ_trans, learning_rate = learning_rates_softmax[iAntenna], activation = cntf.nn.softmax)
        nn_est[:ToepSoftmax] = cntf.ConvNN(nLayers, nAntennas, transform = toep_trans, learning_rate = learning_rates_softmax[iAntenna], activation = cntf.nn.softmax)
    else
        nn_est[:CircReLU]    = cntf.resize(nn_est[:CircReLU],    nAntennas, learning_rate = learning_rates_relu[iAntenna])
        nn_est[:ToepReLU]    = cntf.resize(nn_est[:ToepReLU],    nAntennas, learning_rate = learning_rates_relu[iAntenna])
        nn_est[:CircSoftmax] = cntf.resize(nn_est[:CircSoftmax], nAntennas, learning_rate = learning_rates_softmax[iAntenna])
        nn_est[:ToepSoftmax] = cntf.resize(nn_est[:ToepSoftmax], nAntennas, learning_rate = learning_rates_softmax[iAntenna])
    end

    train!(nn_est, snr = snr, nBatches = nLearningBatches, get_channel = () -> get_channel(nAntennas, nCoherence, nLearningBatchSize), verbose = verbose)

    algs = Dict{Symbol,Any}()
    algs[:GenieMMSE] = (y,h,h_cov) -> mmse_genie(y, h_cov, snr)
    algs[:GenieOMP]  = (y,h,h_cov) -> omp_genie(y, h)
    for (alg,cn) in cn_est
        algs[alg] = (y,h,h_cov) -> mmse.estimate(cn, y)
    end
    for (alg,nn) in nn_est
        algs[alg] = (y,h,h_cov) -> cntf.estimate(nn, y)
    end

    (errs,rates) = evaluate(algs, snr = snr, nBatches = nBatches, get_channel = () -> get_channel(nAntennas, nCoherence, nBatchSize), verbose = verbose)

    for alg in keys(algs)
        new_row = DataFrame(MSE        = errs[alg],
                            rate       = rates[alg],
                            Algorithm  = String(alg),
                            SNR        = snr,
                            nAntennas  = nAntennas,
                            nCoherence = nCoherence)

        if isempty(results)
            results = new_row
        else
            results = vcat(results,new_row)
        end
    end

    if write_file
        CSV.write(filename, results)
    end
    @show results
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
using LazySets, MathematicalPredicates, SparseArrays

function cdplayer_specification()
    # initial set: xᵢ ∈ [-1, 1] for all i
    X0 = BallInf(zeros(120), 1.0)

    # safety property: 2·x₁ - 3·x₂ ≤ 450.8
    property =
        is_contained_in(HalfSpace(sparsevec([1, 2], [2.0, -3.0], 120), 450.8))

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using SparseArrays
using LinearAlgebra
using Random

"""
`hyperkron_graph`
==================
Code to generate a HyperKron Graph using the grass-hopping algorithm.

The original code for this method can be found at:
https://www.cs.purdue.edu/homes/dgleich/codes/hyperkron/

The model definition can be found in the following paper:
https://arxiv.org/abs/1809.03488

Input
-----
- `K`: The initiator tensor
- `k`: The power for the tensor
- `motif`: optional parameter. You can define your own motif to add to the model. Default is undirected triangle.

Functions
---------
-`hyperkron_graph(K,k)`
-`hyperkron_graph(K,k;triangle,False)`

Output
------
- A matrix network type for the hyperkron graph
- A list of hyperedges

Example
--------
A,hedges = hyperkron_graph(kron_params(0.99, 0.2, 0.3, 0.05), 5)
hyperkron_graph(K,k; motif=:triangle, onlytris::Bool=false)

The code calls "hypergraph_to_edges_with_" 'motif' in order
to generate the graph to return. Consequently, you can implement
your own motifs via this routine.
"""
function hyperkron_graph(K,k;
    motif_to_graph::Function = hyperedges_to_graph_with_triangles, kwargs...)

  hedges = fast_hyperkron_edges(K,k; kwargs...)
  return motif_to_graph(size(K,1)^k, hedges), hedges
end

"""
`_next_region!`
===============
Increment to the next region of a Kronecker graph

Input
-----
-`cur`: The current region
-`m`: The total number of elements in the initiated tensor (ie 8 for a 2x2x2 initiator tensor).

Output
------
Outputs the next region
"""
@inline function _next_region!(cur::AbstractArray{T,1}, m::Integer) where T <: Integer
  lastk = length(cur)
  for k=length(cur):-1:1
    lastk = k              # save the last state
    cur[k] += 1            # increment the last element
    if cur[k] == m+1       # if there is spill
      if k > 1             # if there is an array
        continue # recur on prefix
      else
        cur[k] = 0
      end
    end
    break
  end
  for k=lastk+1:length(cur)
    cur[k] = cur[k-1]
  end
  cur
end

"""
`_first_region`
===============
Gives the first region.

Input
-----
-`a`: The vectorized initiator tensor
-`k`: The power of the HyperKron model
Output
------
An array of all ones, representing the first region.
"""
function _first_region(a::AbstractArray,k::Integer)
  return ones(UInt32,k)
end

"""
`num_multiset_permutations`
===========================
Returns the number of multiset permutations of elements in the given array.

Input
-----
-`ms`: A sorted array of numbers, which represent the values in the multiset.
-`start_pos`: Default is 1.

Output
------
The number of multiset permutations.
"""
function num_multiset_permutations(ms::AbstractArray{T,1}, start_pos::Integer = 1) where T<:Integer
  nperm = 1
  count = 1
  countfac = 1
  @inbounds for i in start_pos+1:length(ms)
    nperm *= (i-start_pos+1)
    if ms[i] == ms[i-1]
      count += 1
      countfac *= count
    else
      nperm = div(nperm,countfac)
      # new element
      count = 1
      countfac = 1
    end
  end
  nperm = div(nperm,countfac)
end

"""
`_next_distinct_character`
==========================
Find the next distinct character after a[k]
Returns the index where it occurs or length(a)+1

Input
-----
-`a`: A region in the HyperKron tensor
-`k`: The current character (between 1 and m)
-`m`: Length of a.

Output
------
The index where the next character occurs.
"""
@inline function _next_distinct_character(a, k::Integer, m::Integer)
  # check for easy out
  if k>=m
    return m+1
  else
    ak = a[k]
    if ak != a[k+1] # many cases
      return k+1
    elseif ak == a[m] # some cases
      return m+1
    else
      # we have to search, so there are probably many dups
      # since we assume sorted, we can bisect
      lo = 1
      hi = m-k # we need to search [k+lo:k+hi]
      while hi - lo > 2
        mid = (hi + lo) >> 1
        if ak == a[k+mid]
          lo = mid
        else
          hi = mid
        end
      end
      ak != a[k+lo+1] && return k+lo+1
      return k+lo+2
    end
  end
end

"""
`unrank5!`
==========
This version doesn't use recursion, and unranks the sorted set ms in place
It also gets the multiplicity of each element to update nperm without
recomputing it from scratch, it also can take in the total number
of permutations and just update as we go.
This version tries to avoid integer division by storing nperm as a float

Input
-----
-`ms`: Multiset
-`n`: Current rank
-`nperm`: The number of multiset permutations of the region
-`m`: Length of ms

Output
------
Nothing, it modifies ms.
"""
@inline function unrank5!(ms::AbstractArray{T,1}, n::Integer, nperm, m) where T<:Integer
  multcur = 0
  imultcur = 1.0
  multnext = 0
  npermf = float(nperm)
  @inbounds for k=1:m # for each location in the output
    if n == 1
      break
    end

    if multcur == 0 # update
      multcur = _next_distinct_character(ms, k, m) - k # multiplicity of the current element
      imultcur = 1.0/multcur
    end

    # update nperm for removing multcur
    npermf *= multcur
    npermf /= (m-k+1) # shrink the length
    nperm = round(typeof(nperm),npermf)

    #=
    if nperm != num_multiset_permutations(ms, k+1)
      @show 0, ms, k, num_multiset_permutations(ms, k+1), nperm
    end
    @assert nperm == num_multiset_permutations(ms, k+1)
    =#

    if n <= nperm # this means we start with ms[k]
      multcur -= 1
      if multcur != 0
        imultcur = 1.0/multcur
      end
      continue
    else # we need to search
      n -= nperm

      if k+multcur <= m
        if ms[k] == ms[k+multcur]
          #@show ms, k, _next_distinct_character(ms,k,m), multcur
          @assert ms[k] != ms[k+multcur]
        end
      end

      pos = k+multcur
      while pos <= m
        # we are going to swap ms[k] and ms[pos]
        # but to update, we need to know the multiplicity of ms[pos]

        multnext = _next_distinct_character(ms,pos,m) - pos

        # so the element ms[pos] has multiplicity multnext
        # and element ms[k] has multiplicity multcur-1

        ms[pos], ms[k] = ms[k], ms[pos] # swap leading positions
        npermf *= multnext
        npermf *= imultcur
        nperm = round(typeof(nperm),npermf)

        #=
        if nperm != num_multiset_permutations(ms, k+1)
          @show 1, ms, k, num_multiset_permutations(ms, k+1), nperm
        end
        @assert nperm == num_multiset_permutations(ms, k+1)
        =#

        if n <= nperm
          break
        else
          n -= nperm
          pos += multnext
          multcur = multnext
          imultcur = 1.0/multcur
        end
      end

      multcur = 0 # reset the search
    end
  end
  return ms
end
unrank5!(ms, n) = unrank5!(ms, n, num_multiset_permutations(ms), length(ms))

"""
`_morton_decode3`
=================
Turns an index from a k-dimensional multiplication table into the corresponding index in the 3-dimensional tensor product.

Input
-----
-`mind`: Length k (where k is the number of kronecker products to take) array of indices in the k-dimensional multiplication table
-`n`: The side length of the initiator tensor
-`sz`: The full size of the initiator tensor

Output
------
Length 3 array, representing the corresponding indices in the tensor product.

Example
-------
_morton_decode3([8,8,8,8],2,(2,2,2))
"""
@inline function _morton_decode3(mind::AbstractArray{T,1}, n, sz) where T<:Integer
  row = 0
  col = 0
  tub = 0
  base = 1
  for val in mind
    # CU: ind2sub() was deprecated in Julia 0.7, changed to a new way instead
    s1,s2,s3 = Base._ind2sub(sz, val) #basically backward map where base keeps track of the block size n^
    row += base*(s1-1)
    col += base*(s2-1)
    tub += base*(s3-1)
    base *= n
  end
  return (row+1, col+1, tub+1)
end

#_randgeo(p) = ceil(Int,-randexp() / log1p(-p)) # jsut a copy of Distributions.jl
#_randgeo(p) = ceil(Int,-randexp() / log1p(-p)) # jsut a copy of Distributions.jl
#_randgeo_true(p) = ceil(Int,-randexp() / log1p(-p)) # jsut a copy of Distributions.jl
#_randgeo(p) = _randgeo_true(max(0.01*p, p + 0.0*p*randn()))
_randgeo(p) = ceil(Int,-randexp() / log1p(-p))

"""
`_generate_region!`
==================
Given a region of the tensor, do the grass-hopping procedure to generate hyperedges.

Input
-----
-`ei`: First indices of hyperedges
-`ej`: Second indices of hyperedges
-`ek`: Third indices of hyperedges
-`r`: Region
-`e`: Edge in Region
-`p`: The probability of edges within this particular region.
-`small_n`: The length of one dimension of the initiator matrix.
-`sz`: The full dimensions of the initiator tensor.
-`onlytris`: Boolean
-`onlysym`: Boolean

Output
------
Nothing, modifies input.
"""
function _generate_region!(
    ei::Array{Ta}, ej::Array{Ta}, ek::Array{Ta},
    r::AbstractArray{Tb,1},
    e::AbstractArray{Tb,1},
    p, small_n::Integer, sz,
    onlytris::Bool, onlysym::Bool) where {Ta, Tb <:Integer}

    N = num_multiset_permutations(r) # total area
    m = length(r)
    i = 0

    # ugh, do we need this? should make these tolerances.
    if N*p < 1.0e-12 || p < 1.0e-15
      return
    end

    gap = _randgeo(p)

    while i+gap <= N
      i += gap
      copy!(e, r)
      #@show e, i, N, m
      mult_edge = unrank5!(e, i, N, m) # generate the ith permutation
      #src, dst = _morton_decode(mult_edge, small_n)
      src,dst,tub = _morton_decode3(mult_edge, small_n, sz)
      # only add the edge if it's in the lower region of the matrix.
      # This means we are assuming a symmetric generation.
      if (onlysym && src <= dst <= tub) || !onlysym
        if (onlytris && (src != dst != tub != src)) || !onlytris
          push!(ei, src)
          push!(ej, dst)
          push!(ek, tub)
        end
      end
      gap = _randgeo(p)
    end
end

"""
`fast_hyperkron_edges`
======================
Returns the hyperedges of the Kronecker model

Input
-----
-`K`: Initial tensor for hyperkron model
-`k`: Power for the hyperkron model

Output
------
ei,ej,ek: Three arrays which hold the indices for hyperedges
"""
function fast_hyperkron_edges(K,k;
  noise::Float64=0.0,onlytris::Bool=false,onlysym::Bool=true)
  n = size(K,1)
  sz = size(K)
  v = vec(K) # vectorized

  vn = v/sum(v)
  V = collect(v for _ in 1:k) # create a set of parameters for each level
  #V = collect(max.(0.0,min.(2.0,sum(v).*normalize!(vec(kron_params(vn[1],vn[2]+noise*(2*rand()-1.0),vn[4]+noise*(2*rand()-1.0),vn[8])),1))) for _ in 1:k )
  #@show V

  # Dense = adjacency matrix, every spot is filled in, even if there's a 0. Very large graph = expensive to store
  # Sparse = Store a list of non-zeroes edges. For i and j, there will be two one-dimension array

  # hold edges
  ei = zeros(Int,0)
  ej = similar(ei)
  ek = similar(ei)

  nregions = 0
  # v is the vectoriezd tensor
  region = _first_region(v, k)

  edge_in_region = similar(region) # allocated once to avoid lots of allocs
  while region[1] != 0
    nregions += 1
    # compute the probability of erdos-reyes regions
    p = 1.0
    level = 1
    for j in region
      p *= V[level][j]
      level += 1
    end
    _generate_region!(ei, ej, ek, region, edge_in_region, p, n, sz, onlytris, onlysym)
    _next_region!(region, length(v))
  end
  return ei, ej, ek
end

"""
`kron_params`
=============
Turn 4 parameters into a 2x2x2 symmetric tensor representation for the hyperkron_graph function

Input
-----
-`a`: probability between 0 and 1
-`b`: probability between 0 and 1
-`c`: probability between 0 and 1
-`d`: probability between 0 and 1

Output
Correctly formatted tensor parameters

"""
function kron_params(a,b,c,d)
  K = zeros(2,2,2)
  K[1,1,1] = a
  K[1,1,2] = b; K[1,2,1] = b;   K[2,1,1] = b
  K[1,2,2] = c; K[2,1,2] = c;   K[2,2,1] = c
  K[2,2,2] = d
  return K
end

"""
`hyperedges_to_graph_with_triangles`
====================================
Turns a set of hyperedges into a graph.

Input
-----
-`n`: The size of the graph
-`hedges`: The hyperedges

Output
------
A MatrixNetwork graph based on the given hyperedges.
"""
function hyperedges_to_graph_with_triangles(n::Int, hedges)
  ei = zeros(Int64,0)
  ej = zeros(Int64,0)
  for i=1:length(hedges[1])
    push!(ei, hedges[1][i]); push!(ej, hedges[2][i]) # (i -> j)
    push!(ei, hedges[1][i]); push!(ej, hedges[3][i] ) # (i -> k)
    push!(ei, hedges[2][i]); push!(ej, hedges[3][i] ) # (j -> k)
  end
  A = sparse(ei,ej,1.0,n,n)
  A = max.(A,A') # symmetrize
  A = A - spdiagm(0 => diag(A)) #remove diagonal entries
  fill!(A.nzval, 1.0) # set all non-zeros
  return A
end


text:
for t=1:parameters[:T], (k, v) in results[t][2]
    results[t][2][k] = results[t][2][k][:,:,:,1:1]
end
jldopen("results.jld2", "w") do file
		file["results"] = results
end


text:
#=
- button http://jsxgraph.org/docs/symbols/Button.html
- checkbox http://jsxgraph.org/docs/symbols/Checkbox.html
=#

mutable struct Slider <: Object
    name::String
    vals::Vector{Vector{<:Real}}
    opts::Option{LittleDict{Symbol,Any}}
end

"""
    slider(name, vals; opts...)

Create a new slider, add it to the current board.

* name - name of the slider ex: `"slider1"`
* vals - position and values of the slider ex: `[[0,0],[3,0],[0,1.5,3]]`
          first element is `[xa, ya]` (position of minimum value point)
          second element is `[xb, yb]` (position of maximum value point)
          third element is `[va, v0, vb]` (min value, default value, max value)
"""
function slider(name, vals; kw...)
    @assert length(vals) == 3 "`slider`::expected three subarrays in `vals`."
    @assert length(vals[1]) == 2 && length(vals[2]) == 2 &&
            length(vals[3]) == 3 "`slider`::subvector dims should be 2,2,3."
    s = Slider(name, vals, dict(;kw...))
    return s
end
slider(name, a, b, v; kw...) = slider(name, [a, b, v]; kw...)

# ---------------------------------------------------------------------------

val(s::Slider) = s.vals[3][2] # midvalue of last vector [va, v0, vb]

# ---------------------------------------------------------------------------

function str(s::Slider, b::Board)
    opts = get_opts(s)
    jss = js".create('slider', $(s.vals), $opts);"
    return "$(s.name)=" * b.name * jss.s
end


text:
module LVCyl

export LorentzVectorCyl

export px, py, pz, energy, fast_mass
export deltaphi, deltar, tocartesian, fromcartesian
export ΔR, Δϕ, Δη

import Base: +, *, zero

include("main.jl")

end # module