text:
"""
    filter_vertices_degree_max(vertices, graph::AbstractGraph, degree::Int)

Return the vertices in `vertices` whose degree is lesser than `degree`.
"""
function filter_vertices_degree_max(;vertices, graph::AbstractGraph, degree::Int)
    return filter(vertex -> degreeislesser(graph, vertex, degree), vertices)
end

"""
    filter_vertices_degree_min(vertices, graph::AbstractGraph, degree::Int)

Return the vertices in `vertices` whose degree is greater than `degree`.
"""
function filter_vertices_degree_min(;vertices, graph::AbstractGraph, degree::Int)
    return filter(vertex -> degreeisgreater(graph, vertex, degree), vertices)
end

"""
    filter_vertices_degree(vertices, graph::AbstractGraph, degrees::AbstractVector)

Return the vertices in `vertices` whose degree is in `degrees`.
"""
function filter_vertices_degree(;vertices, graph::AbstractGraph, degrees::AbstractVector)
    return filter(vertex -> degreeis(graph, vertex, degrees), vertices)
end

"""
    filter_vertices_size_biggest_clique_max(vertices, cliques, size::Int)

Return the vertices in `vertices` whose biggest clique size in `cliques` is lesser than `size`.
"""
function filter_vertices_size_biggest_clique_max(;vertices, cliques, size::Int)
    valid_vertices = []
    for vertex in vertices
        filter_cliques = filter(clique -> vertex in clique, cliques)
        if getmaxsize(filter_cliques) <= size
            push!(valid_vertices, vertex)
        end
    end
    return valid_vertices
end

"""
    filter_vertices_size_biggest_clique_min(vertices, cliques, size::Int)

Return the vertices in `vertices` whose biggest clique size in `cliques` is greater than `size`.
"""
function filter_vertices_size_biggest_clique_min(;vertices, cliques, size::Int)
    valid_vertices = []
    for vertex in vertices
        filter_cliques = filter(clique -> vertex in clique, cliques)
        if getmaxsize(filter_cliques) >= size
            push!(valid_vertices, vertex)
        end
    end
    return valid_vertices
end

"""
    filter_vertices_size_biggest_clique(vertices, cliques, sizes::AbstractVector)

Return the vertices in `vertices` whose biggest clique size in `cliques` is `sizes`.
"""
function filter_vertices_size_biggest_clique(;vertices, cliques, sizes::AbstractVector)
    valid_vertices = []
    for vertex in vertices
        filter_cliques = filter(clique -> vertex in clique, cliques)
        if getmaxsize(filter_cliques) in sizes
            push!(valid_vertices, vertex)
        end
    end
    return valid_vertices
end

"""
    filter_vertices_size_smallest_clique_max(vertices, cliques, size::Int)

Return the vertices in `vertices` whose smallest clique size in `cliques` is lesser than `size`.
"""
function filter_vertices_size_smallest_clique_max(;vertices, cliques, size::Int)
    valid_vertices = []
    for vertex in vertices
        filter_cliques = filter(clique -> vertex in clique, cliques)
        if getminsize(filter_cliques)<= size
            push!(valid_vertices, vertex)
        end
    end
    return valid_vertices
end

"""
    filter_vertices_size_smallest_clique_min(;vertices, cliques, size::Int)

Return the vertices in `vertices` whose smallest clique size in `cliques` is greater than `size`.
"""
function filter_vertices_size_smallest_clique_min(;vertices, cliques, size::Int)
    valid_vertices = []
    for vertex in vertices
        filter_cliques = filter(clique -> vertex in clique, cliques)
        if getminsize(filter_cliques)>= size
            push!(valid_vertices, vertex)
        end
    end
    return valid_vertices
end

"""
    filter_vertices_size_smallest_clique(;vertices, cliques, sizes::AbstractVector)

Return the vertices in `vertices` whose smallest clique size in `cliques` is in `sizes`.
"""
function filter_vertices_size_smallest_clique(;vertices, cliques, sizes::AbstractVector)
    valid_vertices = []
    for vertex in vertices
        filter_cliques = filter(clique -> vertex in clique, cliques)
        if getminsize(filter_cliques) in sizes
            push!(valid_vertices, vertex)
        end
    end
    return valid_vertices
end

"""
    filter_vertices_distance_max(vertices, src::Int, graph::AbstractGraph, dist::Int)

Return the vertices in `vertices` whose distance to `src` is lesser than `dist`.
"""
function filter_vertices_distance_max(;vertices, src::Int, graph::AbstractGraph, dist::Int)
    distances = dijkstra_shortest_paths(graph, src).dists
    return filter(vertex -> distances[vertex] <= dist, vertices)
end

"""
    filter_vertices_distance_min(vertices, src::Int, graph::AbstractGraph, dist::Int)

Return the vertices in `vertices` whose distance to `src` is greater than `dist`.
"""
function filter_vertices_distance_min(;vertices, src::Int, graph::AbstractGraph, dist::Int)
    distances = dijkstra_shortest_paths(graph, src).dists
    return filter(vertex -> distances[vertex] >= dist, vertices)
end

"""
    filter_vertices_distance_min(vertices, src::Int, graph::AbstractGraph, dist::AbstractVector)

Return the vertices in `vertices` whose distance to `src` is in `dists`.
"""
function filter_vertices_distance(;vertices, src::Int, graph::AbstractGraph, dists::AbstractVector)
    distances = dijkstra_shortest_paths(graph, src).dists
    return filter(vertex -> distances[vertex] in dists, vertices)
end

"""
    filter_vertices_same_clique(vertices, src::Int, cliques)

Return the vertices in `vertices` that share a similar clique to `src`.
"""
function filter_vertices_same_clique(;vertices, src::Int, cliques)
    return filter(vertex -> !isempty(getsubset(cliques, [src, vertex])), vertices)
end

"""
    filter_vertices_no_same_clique(vertices, src::Int, cliques)

Return the vertices in `vertices` that don't share a similar clique to `src`.
"""
function filter_vertices_no_same_clique(;vertices, src::Int, cliques)
    return filter(vertex -> isempty(getsubset(cliques, [src, vertex])), vertices)
end

"""
    filter_vertices_no_edge(vertices, src::Int, choleskygraph::AbstractGraph)

Return the vertices in `vertices` that don't share an edge with `src`.
"""
function filter_vertices_no_edge(;vertices, src::Int, choleskygraph::AbstractGraph)
    return filter(vertex -> !has_edge(choleskygraph, src, vertex), vertices)
end

"""
    filter_vertices_has_edge(vertices, src::Int, choleskygraph::AbstractGraph)

Return the vertices in `vertices` that share an edge with `src`.
"""
function filter_vertices_has_edge(;vertices, src::Int, choleskygraph::AbstractGraph)
    return filter(vertex -> has_edge(choleskygraph, src, vertex), vertices)
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
module ImageFeatureTests

using ImageFeatures, Images, TestImages, Distributions, ColorTypes
using Test
using LinearAlgebra
import Random.seed!

function check_samples(sample_one, sample_two, size::Int, window::Int)
    check_bool = true
    check_bool = check_bool && length(sample_one) == size
    check_bool = check_bool && length(sample_two) == size
    for s1 in sample_one, s2 in sample_two
        check_bool = check_bool && (s1[1] >= ceil(-window / 2) && s1[1] <= floor((window - 1) / 2)) && (s1[2] >= ceil(-window / 2) && s1[2] <= floor((window - 1) / 2))
        check_bool = check_bool && (s2[1] >= ceil(-window / 2) && s2[1] <= floor((window - 1) / 2)) && (s2[2] >= ceil(-window / 2) && s2[2] <= floor((window - 1) / 2))
    end
    return check_bool
end

function _warp(img, transx, transy)
    res = zeros(eltype(img), size(img))
    for i in 1:size(img, 1) - transx
        for j in 1:size(img, 2) - transy
            res[i + transx, j + transy] = img[i, j]
        end
    end
    res
end

function _warp(img, angle)
    cos_angle = cos(angle)
    sin_angle = sin(angle)
    res = zeros(eltype(img), size(img))
    cx = size(img, 1) / 2
    cy = size(img, 2) / 2
    for i in 1:size(res, 1)
        for j in 1:size(res, 2)
            i_rot = ceil(Int, cos_angle * (i - cx) - sin_angle * (j - cy) + cx)
            j_rot = ceil(Int, sin_angle * (i - cx) + cos_angle * (j - cy) + cy)
            if checkbounds(Bool, img, i_rot, j_rot) res[i, j] = bilinear_interpolation(img, i_rot, j_rot) end
        end
    end
    res
end

function _reverserotate(p, angle, center)
    cos_angle = cos(angle)
    sin_angle = sin(angle)
    return CartesianIndex(floor(Int, sin_angle * (p[2] - center[2]) + cos_angle * (p[1] - center[1]) + center[1]), floor(Int, cos_angle * (p[2] - center[2]) - sin_angle * (p[1] - center[1]) + center[2]))
end

tests = [
    "core.jl",
    "brief.jl",
    "glcm.jl",
    "lbp.jl",
    "corner.jl",
    "orb.jl",
    "freak.jl",
    "brisk.jl",
    "houghtransform.jl",
    "hog.jl"
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end

end

ext: jl


text:
# Rational Quadratic Covariance Function

include("rq_iso.jl")
include("rq_ard.jl")

"""
    RQ(ll::Union{Float64,Vector{Float64}}, lσ::Float64, lα::Float64)

Create Rational Quadratic kernel with length scale `exp.(ll)`, signal standard deviation
`exp(lσ)`, and shape parameter `exp(lα)`.

See also [`RQIso`](@ref) and [`RQArd`](@ref).
"""
RQ(ll::Float64, lσ::Float64, lα::Float64) = RQIso(ll, lσ, lα)
RQ(ll::Vector{Float64}, lσ::Float64, lα::Float64) = RQArd(ll, lσ, lα)


text:
#=
    Helper functions
=#
"""
    get_column(model::LinQuadOptimizer, index::MOI.VariableIndex)

Return the column of the variable `index` in `model`.
"""
function get_column(model::LinQuadOptimizer, index::VarInd)
    return model.variable_mapping[index]
end

"""
    get_column(model::LinQuadOptimizer, variable::MOI.SingleVariable)

Return the column of the variable inside a `MOI.SingleVariable` function in
`model`.
"""
function get_column(model::LinQuadOptimizer, variable::SinVar)
    return get_column(model, variable.variable)
end

#=
    Get variable names
=#
MOI.supports(::LinQuadOptimizer, ::MOI.VariableName, ::Type{VarInd}) = true
MOI.canget(::LinQuadOptimizer, ::MOI.VariableName, ::Type{VarInd}) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.VariableName, index::VarInd)
    if haskey(model.variable_names, index)
        return model.variable_names[index]
    else
        return ""
    end
end

function MOI.set!(model::LinQuadOptimizer, ::MOI.VariableName, index::VarInd, name::String)
    if haskey(model.variable_names_rev, name)
        if model.variable_names_rev[name] != index
            error("Duplicate variable name: $(name)")
        end
    elseif name != ""
        if haskey(model.variable_names, index)
            # we're renaming an existing variable
            old_name = model.variable_names[index]
            delete!(model.variable_names_rev, old_name)
        end
        model.variable_names[index] = name
        model.variable_names_rev[name] = index
    end
    return nothing
end

#=
    Get variable by name
=#
function MOI.canget(model::LinQuadOptimizer, ::Type{MOI.VariableIndex}, name::String)
    return haskey(model.variable_names_rev, name)
end
function MOI.get(model::LinQuadOptimizer, ::Type{MOI.VariableIndex}, name::String)
    return model.variable_names_rev[name]
end

#=
    Get number of variables
=#
MOI.canget(::LinQuadOptimizer, ::MOI.NumberOfVariables) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.NumberOfVariables)
    # TODO(odow): work out why this is not the same as
    # length(model.variable_references)
    return get_number_variables(model)
end

#=
    List of Variable References
=#
MOI.canget(::LinQuadOptimizer, ::MOI.ListOfVariableIndices) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ListOfVariableIndices)
    return model.variable_references
end

#=
    Add a single variable
=#

function MOI.addvariable!(model::LinQuadOptimizer)
    add_variables!(model, 1)
    # assumes we add columns linearly
    model.last_variable_reference += 1
    index = VarInd(model.last_variable_reference)
    model.variable_mapping[index] = MOI.get(model, MOI.NumberOfVariables())
    push!(model.variable_references, index)
    push!(model.variable_primal_solution, NaN)
    push!(model.variable_dual_solution, NaN)
    return index
end

#=
    Add multiple variables
=#

function MOI.addvariables!(model::LinQuadOptimizer, number_to_add::Int)
    previous_vars = MOI.get(model, MOI.NumberOfVariables())
    add_variables!(model, number_to_add)
    variable_indices = VarInd[]
    sizehint!(variable_indices, number_to_add)
    for i in 1:number_to_add
        # assumes we add columns linearly
        model.last_variable_reference += 1
        index = VarInd(model.last_variable_reference)
        push!(variable_indices, index)
        # m.last_variable_reference might not be the column if we have
        # deleted variables.
        model.variable_mapping[index] = previous_vars + i
        push!(model.variable_references, index)
        push!(model.variable_primal_solution, NaN)
        push!(model.variable_dual_solution, NaN)
    end
    return variable_indices
end

#=
    Check if reference is valid
=#

function MOI.isvalid(model::LinQuadOptimizer, index::VarInd)
    if haskey(model.variable_mapping, index)
        column = get_column(model, index)
        if 1 <= column <= MOI.get(model, MOI.NumberOfVariables())
            return true
        end
    end
    return false
end

#=
    Delete a variable
=#

function MOI.delete!(model::LinQuadOptimizer, index::VarInd)
    __assert_valid__(model, index)
    column = get_column(model, index)
    delete_variables!(model, column, column)
    # delete from problem
    deleteat!(model.variable_references, column)
    deleteat!(model.variable_primal_solution, column)
    deleteat!(model.variable_dual_solution, column)
    deleteref!(model.variable_mapping, column, index)
    # delete name if not ""
    if haskey(model.variable_names, index)
        name = model.variable_names[index]
        delete!(model.variable_names_rev, name)
        delete!(model.variable_names, index)
    end
    # Delete any bounds (deleting from a dict without the key does nothing).
    constraint_map = cmap(model)
    delete_from_dictionary_by_value(constraint_map.upper_bound, index)
    delete_from_dictionary_by_value(constraint_map.lower_bound, index)
    delete_from_dictionary_by_value(constraint_map.fixed_bound, index)
    delete_from_dictionary_by_value(constraint_map.interval_bound, index)
    return nothing
end

"""
    delete_from_dictionary_by_value(dict::Dict{K, V}, index::V)

Delete the entries of a dictionary where the value of a `key=>value` pair is
equal to `index`.

This is used to remove bounds when deleting a variable.
"""
function delete_from_dictionary_by_value(dict::Dict{K, V}, index::V) where {K, V}
    for (key, value) in dict
        if value == index
            delete!(dict, key)
        end
    end
end

#=
    MIP starts
=#
function MOI.supports(::LinQuadOptimizer, ::MOI.VariablePrimalStart, ::Type{VarInd})
    return false
end

function MOI.set!(model::LinQuadOptimizer, ::MOI.VariablePrimalStart,
                  indices::Vector{VarInd}, values::Vector{Float64})
    add_mip_starts!(model, get_column.(Ref(model), indices), values)
end

function MOI.set!(model::LinQuadOptimizer, ::MOI.VariablePrimalStart,
                  index::VarInd, value::Float64)
    MOI.set!(model, MOI.VariablePrimalStart(), [index], [value])
end


text:
export Gradient

################################################################################
######   Preconditioning algorithm definition (structure holding params)  ######
################################################################################
"""
    Gradient()

Algorithm for descending along the steepest gradient with SGD-based optimizers.
"""
struct Gradient <: Algorithm end


text:
using Gtk, GLWindow, GLAbstraction, Reactive, GeometryTypes, Colors, GLVisualize
using GtkReactive
using Gtk.GConstants, ModernGL

mutable struct GtkContext <: GLWindow.AbstractContext
    window::Gtk.GLArea
    framebuffer::GLWindow.GLFramebuffer
end
function make_context_current(screen::Screen)
    gl_area = screen.glcontext.window
    Gtk.make_current(gl_area)
end
GLWindow.isopen(x::Gtk.GLArea) = true

global screen = Ref{Screen}()
function init_screen(gl_area, resolution, mesh_color)
    Gtk.make_current(gl_area)
    gtk_context = GtkContext(gl_area, GLWindow.GLFramebuffer(Signal(resolution)))
    window_area = Signal(SimpleRectangle(0, 0, resolution...))
    signals = Dict(
        :mouse_button_released => Reactive.Signal(0),
        :mouse_buttons_pressed => Reactive.Signal(Set(Int[])),
        :scroll => Reactive.Signal(Vec(0.0, 0.0)),
        :buttons_pressed => Reactive.Signal(Set(Int[])),
        :window_size => Reactive.Signal(Vec(resolution...)),
        :window_area => window_area,
        :cursor_position => Reactive.Signal(Vec(0.0,0.0)),
        :mouseinside => Reactive.Signal(true),
        :mouse_button_down => Reactive.Signal(0),
        :mouseposition => Reactive.Signal(Vec(0.0, 0.0)),
        :framebuffer_size => Reactive.Signal(Vec(resolution...)),
        :button_down => Reactive.Signal(0),
        :button_released => Reactive.Signal(0),
        :window_open => Reactive.Signal(true),
        :keyboard_buttons => Reactive.Signal((0,0,0,0)),
        :mouse2id => Signal(GLWindow.SelectionID{Int}(-1, -1))
    )
    screen[] = Screen(
        Symbol("GLVisualize"), window_area, nothing,
        Screen[], signals,
        (), false, true, RGBA(1f0,1f0,1f0,1f0), (0f0, RGBA(0f0,0f0,0f0,0f0)),
        Dict{Symbol, Any}(),
        gtk_context
    )

    GLVisualize.add_screen(screen[])
    timesignal = loop(linspace(0f0,1f0,360))
    rotation_angle  = const_lift(*, timesignal, 2f0*pi)
    start_rotation  = Signal(rotationmatrix_x(deg2rad(90f0)))
    rotation        = map(rotationmatrix_y, rotation_angle)
    final_rotation  = map(*, start_rotation, rotation)
    foreach(final_rotation) do x
        # render a frame each time rotation updates
        Gtk.queue_render(gl_area)
        return
    end
    _view(visualize(loadasset("cat.obj"), color = mesh_color, model = final_rotation), screen[])
    # renderloop(screen[])

end

function render_gtk(window, gtk_area)
    !isopen(window) && return
    fb = GLWindow.framebuffer(window)
    wh = GeometryTypes.widths(window)
    resize!(fb, wh)
    w, h = wh
    #prepare for geometry in need of anti aliasing
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1]) # color framebuffer
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])
    # setup stencil and backgrounds
    glEnable(GL_STENCIL_TEST)
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)
    glStencilMask(0xff)
    glClearStencil(0)
    glClearColor(0,0,0,0)
    glClear(GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT | GL_COLOR_BUFFER_BIT)
    glEnable(GL_SCISSOR_TEST)
    GLWindow.setup_window(window, false)
    glDisable(GL_SCISSOR_TEST)
    glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE)
    # deactivate stencil write
    glEnable(GL_STENCIL_TEST)
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)
    glStencilMask(0x00)
    GLAbstraction.render(window, true)
    glDisable(GL_STENCIL_TEST)

    # transfer color to luma buffer and apply fxaa
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[2]) # luma framebuffer
    glDrawBuffer(GL_COLOR_ATTACHMENT0)
    glClearColor(0,0,0,0)
    glClear(GL_COLOR_BUFFER_BIT)
    glViewport(0, 0, w, h)
    GLAbstraction.render(fb.postprocess[1]) # add luma and preprocess

    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1]) # transfer to non fxaa framebuffer
    glDrawBuffer(GL_COLOR_ATTACHMENT0)
    GLAbstraction.render(fb.postprocess[2]) # copy with fxaa postprocess

    #prepare for non anti aliased pass
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])

    glEnable(GL_STENCIL_TEST)
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)
    glStencilMask(0x00)
    GLAbstraction.render(window, false)
    glDisable(GL_STENCIL_TEST)
    # draw strokes
    glEnable(GL_SCISSOR_TEST)
    GLWindow.setup_window(window, true)
    glDisable(GL_SCISSOR_TEST)
    glViewport(0,0, wh...)
    #Read all the selection queries
    GLWindow.push_selectionqueries!(window)
    Gtk.attach_buffers(gtk_area) # transfer back to window
    glClearColor(0,0,0,0)
    glClear(GL_COLOR_BUFFER_BIT)
    GLAbstraction.render(fb.postprocess[3]) # copy postprocess
    return
end
function connect_renderloop(gl_area, resolution, mesh_color)

end
function setup_screen()
    resolution = (600, 500)
    name = "GTK + GLVisualize"
    parent = Gtk.Window(name, resolution..., true, true)
    Gtk.visible(parent, true)
    Gtk.setproperty!(parent, Symbol("is-focus"), false)
    box = Gtk.Box(:v)
    push!(parent, box)
    mesh_color = Signal(RGBA{Float32}(1,0,0,1))
    gl_area = Gtk.GLArea()
    Gtk.gl_area_set_required_version(gl_area, 3, 3)
    GLAbstraction.new_context()
    push!(box, gl_area)
    Gtk.setproperty!(box, :expand, gl_area, true)
    Gtk.showall(parent)
    Gtk.signal_connect(gl_area, "render") do gl_area, gdk_context
        if !isassigned(screen)
            init_screen(gl_area, resolution, mesh_color)
        end
        render_gtk(screen[], gl_area)
        glFlush()
        return false
    end
    return parent
end

parent = setup_screen();
Gtk.showall(parent)

using Gtk, GLVisualize, GLWindow

function window_area!(window::Screen)
    callback = Gtk.@sigatom function (widget::Gtk.GtkGLArea, width::Int32, height::Int32)
        rect_signal = window.inputs[:window_area]
        push!(rect_signal, IRect(minimum(value(rect_signal)), width, height))
        return true
    end
    gl_area = window.glcontext.window # the Gtk.gl_area
    signal_connect(callback, gl_area, "resize")
    return
end

function scroll_event!(window::Screen)
    callback = Gtk.@sigatom function (widget::Gtk.GtkGLArea, s::Gtk.GdkEventScroll)
        if s.event_type == Gtk.GdkEventType.GDK_SCROLL
            scroll = window.inputs[:scroll]
            dirs = if s.direction == 0
                (0.0, 1.0)
            elseif s.direction == 1
                (0.0, -1.0)
            elseif s.direction == 2
                (-1.0, 0.0)
            elseif s.direction == 3
                (1.0, 0.0)
            else
                (0.0, 0.0)
            end
            push!(scroll, dirs)
            push!(scroll, (0.0, 0.0))
        end
        return true
    end
    gl_area = window.glcontext.window # the Gtk.gl_area
    signal_connect(callback, gl_area, "scroll-event")
    add_events(gl_area, GConstants.GdkEventMask.SCROLL) # same
    return
end

function to_mouse_button(x)
    if x == 1
        0
    elseif x == 2
        1
    elseif x == 3
        2
    else
        # TODO turn into error
        warn("Button is $x, while $(Gtk.GdkModifierType.BUTTON1)")
        Mouse.left
    end
end
function mousebuttons!(window::Screen)
    callback = Gtk.@sigatom function (widget::Gtk.GtkGLArea, event::Gtk.GdkEventButton)
        button = to_mouse_button(event.button)
        action = event.event_type
        if (action in Gtk.GdkEventType.BUTTON_PRESS:Gtk.GdkEventType.TRIPLE_BUTTON_PRESS) ||
                action == Gtk.GdkEventType.BUTTON_RELEASE

            set_s = window.inputs[:mouse_buttons_pressed]
            set = value(set_s)
            if action in Gtk.GdkEventType.BUTTON_PRESS:Gtk.GdkEventType.TRIPLE_BUTTON_PRESS
                push!(set, button)
            elseif action == GdkEventType.BUTTON_RELEASE
                delete!(set, button)
            else
                warn("unknown action: $(action)")
            end
            push!(set_s, set) # trigger event
        end
        return true
    end
    gl_area = window.glcontext.window
    Gtk.add_events(gl_area,
        Gtk.GConstants.GdkEventMask.GDK_BUTTON_PRESS_MASK |
        Gtk.GConstants.GdkEventMask.GDK_BUTTON_RELEASE_MASK
    )
    Gtk.signal_connect(callback, gl_area, "button_press_event")
    Gtk.signal_connect(callback, gl_area, "button_release_event")
    return
end
function mouse!(window::Screen)
    callback = Gtk.@sigatom function (widget::Gtk.GtkGLArea, s::Gtk.GdkEventMotion)
        if s.event_type == Gtk.GdkEventType.GDK_MOTION_NOTIFY
            pos = window.inputs[:mouseposition]
            push!(pos, (s.x, s.y))
        end
        return true
    end
    gl_area = window.glcontext.window # the Gtk.gl_area
    add_events(gl_area, GConstants.GdkEventMask.POINTER_MOTION)
    signal_connect(callback, gl_area, "motion-notify-event")
    return true
end
scroll_event!(screen[])
mouse!(screen[])
mousebuttons!(screen[])


text:
using Test, BellScenario

@testset "./src/verifications.jl" begin

using SignalingDimension

@testset "verify_facet()" begin
    @testset "3-2-3 signaling polytope examples" begin
        facet = BellGame([1 0 0;0 1 0;0 0 1], 2)
        aff_ind_vertices = [
            [1 0 1;0 1 0;0 0 0], [1 1 0;0 0 0;0 0 1],
            [0 0 0;1 1 0;0 0 1], [1 0 0;0 1 1;0 0 0],
            [1 0 0;0 0 0;0 1 1], [0 0 0;0 1 0;1 0 1]
        ]

        @test verify_facet(2, facet, aff_ind_vertices)

        @test !verify_facet(1, facet, aff_ind_vertices)
        @test !verify_facet(3, facet, aff_ind_vertices)

        @test !verify_facet(2, facet, aff_ind_vertices[1:5])
        @test !verify_facet(2, facet, [aff_ind_vertices...,[0 0 0;0 1 0;1 0 1]])
    end

    @testset "invalid vertices examples" begin
        facet = BellGame([1 0 0;0 1 0;0 0 0], 2)
        aff_ind_vertices = [
            [1 0 1;0 1 0;0 0 0], [1 1 0;0 1 0;0 0 1],
            [1 0 0;1 1 0;0 0 1], [1 0 0;0 1 0;0 0 0],
            [1 0 0;0 1 0;0 1 1], [1 0 0;0 1 0;1 0 1]
        ]

        @test LocalPolytope.dimension(aff_ind_vertices) == 5
        @test !verify_facet(2, facet, aff_ind_vertices)
    end
end

end


text:
using Distributions
using LinearAlgebra
using Random
using Plots
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.ParameterDistributionStorage

# Seed for pseudo-random number generator for reproducibility
rng_seed = 41
Random.seed!(rng_seed)

# Number of synthetic observations from G(u)
n_obs = 1
# Defining the observation noise level
noise_level =  1e-8
# Independent noise for synthetic observations
Γy = noise_level * Matrix(I, n_obs, n_obs)
noise = MvNormal(zeros(n_obs), Γy)

# Loss Function (unique minimum)
function G(u)
    return [sqrt((u[1]-1)^2 + (u[2]+1)^2)]
end

# Loss Function Minimum
u_star = [1.0, -1.0]
y_obs  = G(u_star) + 0 * rand(noise)

# Define Prior
prior_distns = [Parameterized(Normal(0., sqrt(1))),
                Parameterized(Normal(-0., sqrt(1)))]
constraints = [[no_constraint()], [no_constraint()]]
prior_names = ["u1", "u2"]
prior = ParameterDistribution(prior_distns, constraints, prior_names)
prior_mean = get_mean(prior)
prior_cov = get_cov(prior)

# Calibrate
N_ens = 50  # number of ensemble members
N_iter = 20 # number of EKI iterations
initial_ensemble = EnsembleKalmanProcessModule.construct_initial_ensemble(prior, N_ens;
                                                rng_seed=rng_seed)

ekiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_ensemble,
                    y_obs, Γy, Inversion())
#
for i in 1:N_iter
    params_i = get_u_final(ekiobj)
    g_ens = hcat([G(params_i[:,i]) for i in 1:N_ens]...)
    EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens)
end

u_init = get_u_prior(ekiobj)

for i in 1:N_iter
    u_i = get_u(ekiobj,i)
    p = plot(u_i[1,:], u_i[2,:], seriestype=:scatter, xlims = extrema(u_init[1,:]), ylims = extrema(u_init[2,:]))
    plot!([u_star[1]], xaxis="u1", yaxis="u2", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "EKI iteration = " * string(i)
        )
    plot!([u_star[2]], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optimum")
    display(p)
    sleep(0.1)
end

##
rng_seed = 10 # 10 converges to one minima 100 converges to the other

# Loss Function (two minima)
function G(u)
    return [abs((u[1]-1)*(u[1]+1))^2 + (u[2]+1)^2]
end

# Loss Function Minimum
u_star1 = [1.0, -1.0]
u_star2 = [-1.0, -1.0]
G(u_star1)[1] == G(u_star2)[1]
y_obs  = [0.0]

# Define Prior
prior_distns = [Parameterized(Normal(0., sqrt(2))),
                Parameterized(Normal(-0., sqrt(2)))]
constraints = [[no_constraint()], [no_constraint()]]
prior_names = ["u1", "u2"]
prior = ParameterDistribution(prior_distns, constraints, prior_names)
prior_mean = get_mean(prior)
prior_cov = get_cov(prior)

# Calibrate
N_ens = 50  # number of ensemble members
N_iter = 40 # number of EKI iterations
initial_ensemble = EnsembleKalmanProcessModule.construct_initial_ensemble(prior, N_ens;
                                                rng_seed=rng_seed)

ekiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_ensemble,
                    y_obs, Γy, Inversion())
#
for i in 1:N_iter
    params_i = get_u_final(ekiobj)
    g_ens = hcat([G(params_i[:,i]) for i in 1:N_ens]...)
    EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens)
end

u_init = get_u_prior(ekiobj)
for i in 1:N_iter
    u_i = get_u(ekiobj,i)
    p = plot(u_i[1,:], u_i[2,:], seriestype=:scatter, xlims = (-2,2), ylims = (-2,2))
    plot!([1], xaxis="u1", yaxis="u2", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "EKI iteration = " * string(i)
        )
    plot!([-1], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optima 1")

    plot!([-1], xaxis="u1", yaxis="u2", seriestype="vline",
        linestyle=:dash, linecolor=:green, label = false,
        title = "EKI iteration = " * string(i)
        )
    plot!([-1], seriestype="hline", linestyle=:dash, linecolor=:green, label = "optima 2")
    display(p)
    sleep(0.1)
end


text:
using Combinatorics: permutations

export non0hist, genentropy, permentropy

"""
```julia
non0hist(ε, dataset::AbstractDataset)
```
Partition a dataset into tabulated intervals (boxes) of
size `ε` and return the sum-normalized histogram in an unordered 1D form,
discarding all zero elements.

## Performances Notes
This method is effecient in both memory
and speed, because it uses a dictionary to collect the information of bins with
elements, while it completely disregards empty bins. This allows
computation of entropies of high-dimensional datasets and
with small box sizes `ε` without memory overflow.

Use e.g. `fit(Histogram, ...)` from
[`StatsBase`](http://juliastats.github.io/StatsBase.jl/stable/) if you
wish to keep information about the edges of the binning as well
as the zero elements.
"""
function non0hist(ε::Real, data::AbstractDataset{D, T}) where {D, T<:Real}
    # Initialize:
    mini = minima(data)
    L = length(data)
    # Perform "Histogram":
    # `d` is a dictionary that contains all the histogram information
    # the keys are the bin edges indices and the values are the amount of
    # datapoints in each bin
    d = Dict{SVector{D, Int}, Int}()
    for point in data
        # index of datapoint in the ranges space:
        # It is necessary to convert floor to Int (representation issues)
        ind::SVector{D, Int} = floor.(Int, (point - mini)/ε)

        # Add 1 to the bin that contains the datapoint:
        d[ind] = get(d, ind, 0) + 1
    end
    return collect(values(d))/L
end

"""
    genentropy(α, ε, dataset::AbstractDataset; base = e)
Compute the `α` order generalized (Rényi) entropy [1] of a dataset,
by first partitioning it into boxes of length `ε` using [`non0hist`](@ref).
```julia
genentropy(α, p::AbstractArray; base = e)
```
Compute the entropy of an array `p` directly, assuming that `p` is
sum-normalized.

Optionally use `base` for the logarithms.

## Description
Let ``p`` be an array of probabilities (summing to 1). Then the Rényi entropy is
```math
R_\\alpha(p) = \\frac{1}{1-\\alpha} \\log \\left(\\sum_i p[i]^\\alpha\\right)
```
and generalizes other known entropies,
like e.g. the information entropy
(``\\alpha = 1``, see [2]), the maximum entropy (``\\alpha=0``,
also known as Hartley entropy), or the correlation entropy
(``\\alpha = 2``, also known as collision entropy).

## References

[1] : A. Rényi, *Proceedings of the fourth Berkeley Symposium on Mathematics,
Statistics and Probability*, pp 547 (1960)

[2] : C. E. Shannon, Bell Systems Technical Journal **27**, pp 379 (1948)
"""
function genentropy(α::Real, ε::Real, data::AbstractDataset;
    base=Base.MathConstants.e)
    ε ≤ 0 && throw(ArgumentError("Box-size for entropy calculation must be > 0."))
    p = non0hist(ε, data)
    return genentropy(α, p; base = base)
end
genentropy(α::Real, ε::Real, matrix; base = Base.MathConstants.e) =
genentropy(α, ε, Dataset(matrix); base = base)

function genentropy(α::Real, p::AbstractArray{T}; base = e) where {T<:Real}
  α < 0 && throw(ArgumentError("Order of Rényi entropy must be ≥ 0."))

  if α ≈ 0
    return log(base, length(p)) #Hartley entropy, max-entropy
  elseif α ≈ 1
    return -sum( x*log(base, x) for x in p ) #Shannon entropy
  elseif isinf(α)
    return -log(base, maximum(p)) #Min entropy
  else
    return (1/(1-α))*log(base, sum(x^α for x in p) ) #Renyi α entropy
  end
end

"""
    permentropy(x::AbstractVector, order [, interval=1]; base = e)

Compute the permutation entropy [1] of given `order`
from the `x` timeseries.

Optionally, `interval` can be specified to
use `x[t0:interval:t1]` when calculating permutation of the
sliding windows between `t0` and `t1 = t0 + interval * (order - 1)`.

Optionally use `base` for the logarithms.

## References

[1] : C. Bandt, & B. Pompe, [Phys. Rev. Lett. **88** (17), pp 174102 (2002)](http://doi.org/10.1103/PhysRevLett.88.174102)
"""
function permentropy(
        time_series::AbstractArray{T, 1}, orderi::Integer,
        interval::Integer = 1;
        base=Base.MathConstants.e) where {T}

    orderi > 255 && throw(ArgumentError("order = $orderi is too large, "*
              "must be smaller than $(Int(typemax(UInt8))) can be used."))
    order = UInt8(orderi)
    # To use `searchsortedfirst`, we need each permutation to be
    # "comparable" (ordered) type.  Let's use NTuple here:
    PermType = NTuple{Int(order), UInt8}

    perms = map(PermType, permutations(1:order))
    count = zeros(UInt64, length(perms))

    for t in 1:length(time_series) - interval * order + 1
        sample = @view time_series[t:interval:t + interval * (order - 1)]
        i = searchsortedfirst(perms, PermType(sortperm(sample)))
        count[i] += 1
    end

    # To compute `p log(p)` correctly for `p = 0`, we first discard
    # cases with zero occurrence.  They don't contribute to the final
    # sum hence to the entropy:
    nonzero = [c for c in count if c != 0]

    p = nonzero ./ sum(nonzero)
    return -sum(p .* log.(base, p))
end


text:
"""
    SimpleWeightedDiGraph{T, U}

A type representing a directed graph with weights of type `U`.
"""
mutable struct SimpleWeightedDiGraph{T<:Integer, U<:Real} <: AbstractSimpleWeightedGraph{T, U}
    weights::SparseMatrixCSC{U, T} # indexed by [dst, src]
end

ne(g::SimpleWeightedDiGraph) = nnz(g.weights)

# Graph{UInt8}(6), Graph{Int16}(7), Graph{UInt8}()
function (::Type{SimpleWeightedDiGraph{T, U}})(n::Integer = 0) where T<:Integer where U<:Real
    weights = spzeros(U, T, T(n), T(n))
    return SimpleWeightedDiGraph{T, U}(weights)
end

# Graph()
SimpleWeightedDiGraph() = SimpleWeightedDiGraph{Int, Float64}()

# Graph(6), Graph(0x5)
SimpleWeightedDiGraph(n::T) where T<:Integer = SimpleWeightedDiGraph{T, Float64}(n)

# Graph(UInt8)
SimpleWeightedDiGraph(::Type{T}) where T<:Integer = SimpleWeightedDiGraph{T, Float64}(zero(T))

# Graph(UInt8, Float32)
SimpleWeightedDiGraph(::Type{T}, ::Type{U}) where T<:Integer where U<:Real = SimpleWeightedDiGraph{U, T}(zero(T))

# DiGraph(AbstractSimpleGraph)
function SimpleWeightedDiGraph(g::LightGraphs.SimpleGraphs.AbstractSimpleGraph, ::Type{U}=Float64) where U <: Real
    T = eltype(g)
    return SimpleWeightedDiGraph{T, U}(adjacency_matrix(g)')
end

# DiGraph(AbstractSimpleGraph, defaultweight)
function SimpleWeightedDiGraph(g::LightGraphs.SimpleGraphs.AbstractSimpleGraph, x::U) where U <: Real
    T = eltype(g)
    return SimpleWeightedDiGraph{T, U}(x.*adjacency_matrix(g, U)')
end

# SimpleWeightedGraph{T, U}(SimpleGraph)
function (::Type{SimpleWeightedDiGraph{T, U}})(g::LightGraphs.SimpleGraphs.SimpleDiGraph)  where T<:Integer where U <: Real
    SimpleWeightedDiGraph{T, U}(adjacency_matrix(LightGraphs.SimpleGraphs.SimpleDiGraph{T}(g), U))
end

# DiGraph(srcs, dsts, weights)
function SimpleWeightedDiGraph(i::AbstractVector{T}, j::AbstractVector{T}, v::AbstractVector{U}; combine = +) where T<:Integer where U<:Real
    m = max(maximum(j), maximum(i))
    SimpleWeightedDiGraph{T, U}(sparse(j, i, v, m, m, combine))
end

# Graph{UInt8}(adjmx)
# function (::Type{SimpleWeightedDiGraph{T, U}})(adjmx::AbstractMatrix) where T<:Integer where U <: Real
#     dima,dimb = size(adjmx)
#     isequal(dima,dimb) || error("Adjacency / distance matrices must be square")
#     issymmetric(adjmx) || error("Adjacency / distance matrices must be symmetric")
#     g = SimpleWeightedDiGraph(U.(LinearAlgebra.fillstored!(copy(adjmx), 1)))
# end

# converts Graph{Int} to Graph{Int32}
# function (::Type{SimpleWeightedDiGraph{T, U}})(g::SimpleWeightedDiGraph) where T<:Integer where U<:Real
#     h_fadj = [Vector{T}(x) for x in fadj(g)]
#     return SimpleGraph(ne(g), h_fadj)
# end

# Graph(adjmx)
SimpleWeightedDiGraph(adjmx::AbstractMatrix) = SimpleWeightedDiGraph{Int, eltype(adjmx)}(adjmx')

# Graph(digraph). Weights will be added.
# TODO: uncomment this.
SimpleWeightedDiGraph(g::AbstractSimpleWeightedGraph) = SimpleWeightedDiGraph(g.weights)

edgetype(::SimpleWeightedDiGraph{T, U}) where T<:Integer where U<:Real = SimpleWeightedGraphEdge{T,U}

edges(g::SimpleWeightedDiGraph) = (SimpleWeightedEdge(x[2], x[1], x[3]) for x in zip(findnz(g.weights)...))
weights(g::SimpleWeightedDiGraph) = g.weights'
function inneighbors(g::SimpleWeightedDiGraph)
    mat = SparseMatrixCSC(g.weights')
    return [mat.rowval[mat.colptr[i]:mat.colptr[i + 1] - 1] for i in 1:nv(g)]
end

function inneighbors(g::SimpleWeightedDiGraph, v::Integer)
    mat = SparseMatrixCSC(g.weights')
    return mat.rowval[mat.colptr[v]:mat.colptr[v + 1] - 1]
end

# add_edge! will overwrite weights.
function add_edge!(g::SimpleWeightedDiGraph, e::SimpleWeightedGraphEdge)
    @warn "Note: adding edges to this graph type is not performant." maxlog=1 _id=:swd_add_edge
    T = eltype(g)
    U = weighttype(g)
    s_, d_, w = Tuple(e)

    if w == zero(U)
        @warn "Note: adding edges with a zero weight to this graph type has no effect." maxlog=1 _id=:swd_add_edge_zero
        return false
    end

    s = T(s_)
    d = T(d_)
    (s in vertices(g) && d in vertices(g)) || return false
    g.weights[d, s] = w
    return true
end

function rem_edge!(g::SimpleWeightedDiGraph, e::SimpleWeightedGraphEdge)
    has_edge(g, e) || return false
    U = weighttype(g)
    g.weights[dst(e), src(e)] = zero(U)
    return true
end

copy(g::SimpleWeightedDiGraph) =  SimpleWeightedDiGraph(copy(g.weights))

==(g::SimpleWeightedDiGraph, h::SimpleWeightedDiGraph) = g.weights == h.weights

"""
    is_directed(g)

Return `true` if `g` is a directed graph.
"""
is_directed(::Type{SimpleWeightedDiGraph}) = true
is_directed(::Type{SimpleWeightedDiGraph{T,U}}) where T where U = true
is_directed(g::SimpleWeightedDiGraph) = true


text:
"""
bicgstab_test_pde(n)

PDE test from FR16. Test of kl_bicgstab with all kinds of preconditioning.
"""
function bicgstab_test_pde(n; write = false, eta = 9.8 * 1.e-4)
    pdata = pdegminit(n)
    RHS = pdata.RHS
    ue = pdata.ue
    u0 = zeros(n * n)
#
    # Solve with left preconditioning hard-wired in
#
    fdata = pdata.fdata
    b = Pfish2d(RHS, fdata)
    goutp = kl_bicgstab(u0, b, pdelpatv, eta; pdata = pdata, lmaxit=200)
    pcres = goutp.reshist
    pcres /= pcres[1]
    sollhw = goutp.sol
#
    # Solve with right preconditioning hard-wired in
#
    goutrp = kl_bicgstab(u0, RHS, pderatv, eta; pdata = pdata, lmaxit=200)
    pcresr = goutrp.reshist
    pcresr /= pcresr[1]
    solrhw = copy(u0)
    solrhw .= pdeptv(goutrp.sol,pdata)
    solldiff=norm(solrhw-sollhw,Inf)
#
    # Solve with right preconditioning
#
    goutrp1 = kl_bicgstab(u0, RHS, pdeatv, eta, pdeptv;
             pdata = pdata, lmaxit=200)
    pcresr1 = goutrp1.reshist
    pcresr1 /= pcresr1[1]
    solr = goutrp1.sol
#    solldiff += norm(solr-sollhw,Inf)
    solldiff = max(solldiff,norm(solr-sollhw,Inf))
#
    # Solve with left preconditioning
#
    goutl1 = kl_bicgstab(u0, RHS, pdeatv, eta, pdeptv;
             pdata = pdata, lmaxit=200, side="left")
    pcresl1 = goutl1.reshist
    pcresl1 /= pcresl1[1]
    soll = goutl1.sol
#    solldiff += norm(soll-sollhw,Inf)
    solldiff = max(solldiff,norm(soll-sollhw,Inf))
#
# Hardwired and normal give same results?
#
    leftdel= norm(soll-sollhw,Inf) + norm(pcres-pcresl1,Inf)
    leftpass=(leftdel < 1.e-15)
    rightdel = norm(solr-solrhw,Inf) + norm(pcresr1-pcresr,Inf)
    rightpass=(rightdel < 1.e-15)
#
    # Solve with no preconditioning to duplicate fig 3.4 in red book
#
    goutnp = kl_bicgstab(u0, RHS, pdeatv, eta; lmaxit=200, pdata = pdata)
    pcresnp = goutnp.reshist
    pcresnp /= pcresnp[1]
    solnone=goutnp.sol
    solldiff = max(solldiff,norm(solnone-sollhw,Inf))
#
# Solve with no preconditioning to get a failure
#
    goutnf = kl_bicgstab(u0, RHS, pdeatv, V, eta; lmaxit=20, pdata = pdata)
    failpass = ~goutnf.idid && (goutnf.lits==20)
#
# Are the answers close enough?
#
    println(solldiff/eta)
    sollpass=(solldiff < 2.0*eta)
#
    pass=sollpass && rightpass && leftpass && failpass
    return pass
end


text:
using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test,
        MatrixFactorizations, ArrayLayouts, LinearAlgebra, Random, LazyBandedMatrices, StaticArrays
import InfiniteLinearAlgebra: qltail, toeptail, tailiterate , tailiterate!, tail_de, ql_X!,
                    InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, InfBandCartesianIndices,
                    rightasymptotics, QLHessenberg, ConstRows, PertConstRows, chop, chop!,
                    BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout
import Base: BroadcastStyle
import BlockArrays: _BlockArray
import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
import MatrixFactorizations: QLPackedQ
import BandedMatrices: bandeddata, _BandedMatrix, BandedStyle
import LazyArrays: colsupport, MemoryLayout, ApplyLayout, LazyArrayStyle, arguments
import InfiniteArrays: OneToInf, oneto, RealInfinity
import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayout, LazyBandedLayout

@testset "chop" begin
    a = randn(5)
    b = [a; zeros(5)]
    chop!(b, eps())
    @test b == a

    @test isempty(chop!([0]))

    A = randn(5,5)
    @test chop([A zeros(5,2); zeros(2,5) zeros(2,2)],eps()) == A

    c = PseudoBlockArray([randn(5); zeros(10)], (blockedrange(1:5),))
    d = chop!(c, 0);
    @test length(d) == 6
end

include("test_infconv.jl")
include("test_infbanded.jl")

@testset "∞-block arrays" begin
    @testset "fixed block size" begin
        k = Base.OneTo.(oneto(∞))
        n = Fill.(oneto(∞),oneto(∞))
        @test broadcast(length,k) ≡ map(length,k) ≡ OneToInf()
        @test broadcast(length,n) ≡ map(length,n) ≡ OneToInf()

        b = mortar(Fill([1,2],∞))
        @test blockaxes(b,1) ≡ Block.(OneToInf())
        @test b[Block(5)] == [1,2]
        @test b[Block.(2:∞)][Block.(2:10)] == b[Block.(3:11)]
        @test exp.(b)[Block.(2:∞)][Block.(2:10)] == exp.(b[Block.(3:11)])

        c = PseudoBlockArray(1:∞,Vcat(2,Fill(3,∞)))
        @test c[Block.(2:∞)][Block.(2:10)] == c[Block.(3:11)]

        @test length(axes(b,1)) ≡ ℵ₀
        @test last(axes(b,1)) ≡ RealInfinity()
        @test Base.BroadcastStyle(typeof(b)) isa LazyArrayStyle{1}

        @test unitblocks(oneto(∞)) ≡ blockedrange(Ones{Int}(∞))
        @test unitblocks(2:∞) == 2:∞

        @test unitblocks(oneto(∞))[Block.(2:∞)] == 2:∞
    end

    @testset "1:∞ blocks" begin
        a = blockedrange(oneto(∞))
        @test axes(a,1) == a
        o = Ones((a,))
        @test Base.BroadcastStyle(typeof(a)) isa LazyArrayStyle{1}
        b = exp.(a)
        @test axes(b,1) == a
        @test o .* b isa typeof(b)
        @test b .* o isa typeof(b)
    end

    @testset "concat" begin
        a = unitblocks(1:∞)
        b = exp.(a)
        c = BlockBroadcastArray(vcat,a,b)
        @test length(c) == ∞
        @test blocksize(c) == (∞,)
        @test c[Block(5)] == [a[5],b[5]]

        A = unitblocks(BandedMatrix(0 => 1:∞, 1=> Fill(2.0,∞), -1 => Fill(3.0,∞)))
        B = BlockBroadcastArray(hvcat, 2, A, Zeros(axes(A)), Zeros(axes(A)), A)
        @test B[Block(3,3)] == [A[3,3] 0; 0 A[3,3]]
        @test B[Block(3,4)] == [A[3,4] 0; 0 A[3,4]]
        @test B[Block(3,5)] == [A[3,5] 0; 0 A[3,5]]
        @test blockbandwidths(B) == (1,1)
        @test subblockbandwidths(B) == (0,0)
        @test B[Block.(1:10), Block.(1:10)] isa BlockSkylineMatrix

        C = BlockBroadcastArray(hvcat, 2, A, A, A, A)
        @test C[Block(3,3)] == fill(A[3,3],2,2)
        @test C[Block(3,4)] == fill(A[3,4],2,2)
        @test C[Block(3,5)] == fill(A[3,5],2,2)
        @test blockbandwidths(C) == (1,1)
        @test subblockbandwidths(C) == (1,1)
        @test B[Block.(1:10), Block.(1:10)] isa BlockSkylineMatrix
    end

    @testset "KronTrav" begin
        Δ = BandedMatrix(1 => Ones(∞), -1 => Ones(∞))/2
        A = KronTrav(Δ - 2I, Eye(∞))
        @test axes(A,1) isa InfiniteLinearAlgebra.OneToInfBlocks
        V = view(A, Block.(Base.OneTo(3)), Block.(Base.OneTo(3)))
        @test MemoryLayout(V) isa BlockBandedMatrices.BandedBlockBandedLayout

        u = A * [1; zeros(∞)]
        @test u[1:3] == A[1:3,1]
        @test bandwidths(view(A, Block(1,1))) == (1,1)
    end

    @testset "triangle recurrences" begin
        @testset "n and k" begin
            n = mortar(Fill.(oneto(∞),oneto(∞)));
            k = mortar(Base.OneTo.(oneto(∞)));

            @test n[Block(5)] ≡ layout_getindex(n, Block(5)) ≡ view(n,Block(5)) ≡ Fill(5,5)
            @test k[Block(5)] ≡ layout_getindex(k, Block(5)) ≡ view(k,Block(5)) ≡ Base.OneTo(5)
            @test Base.BroadcastStyle(typeof(n)) isa LazyArrays.LazyArrayStyle{1}
            @test Base.BroadcastStyle(typeof(k)) isa LazyArrays.LazyArrayStyle{1}

            N = 1000
            v = view(n,Block.(Base.OneTo(N)))
            @test view(v,Block(2)) ≡ Fill(2,2)
            @test axes(v) isa Tuple{BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,Base.OneTo{Int64}}}}
            @test @allocated(axes(v)) ≤ 40

            dest = PseudoBlockArray{Float64}(undef, axes(v))
            @test copyto!(dest, v) == v
            @test @allocated(copyto!(dest, v)) ≤ 40

            v = view(k,Block.(Base.OneTo(N)))
            @test view(v,Block(2)) ≡ Base.OneTo(2)
            @test axes(v) isa Tuple{BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,Base.OneTo{Int64}}}}
            @test @allocated(axes(v)) ≤ 40
            @test copyto!(dest, v) == v

            @testset "stack overflow" begin
                i = Base.to_indices(k, (Block.(2:∞),))[1].indices;
                last(i)
            end

            v = view(k,Block.(2:∞))
            @test Base.BroadcastStyle(typeof(v)) isa LazyArrayStyle{1}
            @test v[Block(1)] == 1:2
            @test v[Block(1)] ≡ k[Block(2)] ≡ Base.OneTo(2)

            @test axes(n,1) isa BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,OneToInf{Int64}}}
        end

        @testset "BlockHcat copyto!" begin
            n = mortar(Fill.(oneto(∞),oneto(∞)))
            k = mortar(Base.OneTo.(oneto(∞)))

            a = b = c = 0.0
            dat = BlockHcat(
                BroadcastArray((n,k,b,bc1) -> (k + b-1) * (n + k + bc1) / (2k + bc1), n, k, b, b+c-1),
                BroadcastArray((n,k,abc,bc,bc1) -> (n + k + abc) * (k + bc) / (2k + bc1), n, k, a+b+c,b+c,b+c-1)
                )
            N = 1000
            KR = Block.(Base.OneTo(N))
            V = view(dat,Block.(Base.OneTo(N)),:)
            @test MemoryLayout(V) isa LazyArrays.ApplyLayout{typeof(hcat)}
            @test PseudoBlockArray(V)[Block.(1:5),:] == dat[Block.(1:5),:]
            V = view(dat',:,Block.(Base.OneTo(N)))
            @test MemoryLayout(V) isa LazyArrays.ApplyLayout{typeof(vcat)}
            a = dat.arrays[1]'
            N = 100
            KR = Block.(Base.OneTo(N))
            v = view(a,:,KR)
            @time r = PseudoBlockArray(v)
            @test v == r
        end

        @testset "BlockBanded" begin
            a = b = c = 0.0
            n = mortar(Fill.(oneto(∞),oneto(∞)))
            k = mortar(Base.OneTo.(oneto(∞)))
            Dy = BlockBandedMatrices._BandedBlockBandedMatrix((k .+ (b+c))', axes(k,1), (-1,1), (-1,1))
            N = 100;
            @test Dy[Block.(1:N), Block.(1:N)] == BlockBandedMatrices._BandedBlockBandedMatrix((k .+ (b+c))[Block.(1:N)]', axes(k,1)[Block.(1:N)], (-1,1), (-1,1))
        end
    end

    @testset "blockdiag" begin
        D = Diagonal(mortar(Fill.((-(0:∞)-(0:∞).^2), 1:2:∞)))
        x = [randn(5); zeros(∞)]
        x̃ = PseudoBlockArray(x, (axes(D,1),))
        @test (D * x)[1:10] == (D * x̃)[1:10]
    end

    @testset "sortedunion" begin
        a = cumsum(1:2:∞)
        @test BlockArrays.sortedunion(a,a) ≡ a
        @test BlockArrays.sortedunion([∞],a) ≡ BlockArrays.sortedunion(a,[∞]) ≡ a
        @test BlockArrays.sortedunion([∞],[∞]) == [∞]

        b = Vcat([1,2],3:∞)
        c = Vcat(1,3:∞)
        @test BlockArrays.sortedunion(b,b) ≡ b
        @test BlockArrays.sortedunion(c,c) ≡ c
    end
end

@testset "Algebra" begin
    @testset "BandedMatrix" begin
        A = BandedMatrix(-3 => Fill(7/10,∞), -2 => 1:∞, 1 => Fill(2im,∞))
        @test A isa BandedMatrix{ComplexF64}
        @test A[1:10,1:10] == diagm(-3 => Fill(7/10,7), -2 => 1:8, 1 => Fill(2im,9))

        A = BandedMatrix(0 => Vcat([1,2,3],Zeros(∞)), 1 => Vcat(1, Zeros(∞)))
        @test A[1,2] == 1

        A = BandedMatrix(-3 => Fill(7/10,∞), -2 => Fill(1,∞), 1 => Fill(2im,∞))
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'
        @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])

        A = BandedMatrix(-1 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
        @test MemoryLayout(typeof(view(A.data,:,1:10))) == ApplyLayout{typeof(hcat)}()
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'
        @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])

        A = BandedMatrix(-2 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im,2,3],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'
        @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])

        A = _BandedMatrix(Fill(1,4,∞),ℵ₀,1,2)
        @test A^2 isa BandedMatrix
        @test (A^2)[1:10,1:10] == (A*A)[1:10,1:10] == (A[1:100,1:100]^2)[1:10,1:10]
        @test A^3 isa ApplyMatrix{<:Any,typeof(*)}
        @test (A^3)[1:10,1:10] == (A*A*A)[1:10,1:10]  == ((A*A)*A)[1:10,1:10]  == (A*(A*A))[1:10,1:10] == (A[1:100,1:100]^3)[1:10,1:10]

        @testset "∞ x finite" begin
            A = BandedMatrix(1 => 1:∞) + BandedMatrix(-1 => Fill(2,∞))
            B = _BandedMatrix(randn(3,5), ℵ₀, 1,1)

            @test lmul!(2.0,copy(B)')[:,1:10] ==  (2B')[:,1:10]

            @test_throws ArgumentError BandedMatrix(A)
            @test A*B isa MulMatrix
            @test B'A isa MulMatrix

            @test all(diag(A[1:6,1:6]) .=== zeros(6))

            @test (A*B)[1:7,1:5] ≈ A[1:7,1:6] * B[1:6,1:5]
            @test (B'A)[1:5,1:7] ≈ (B')[1:5,1:6] * A[1:6,1:7]
        end
    end

    @testset "BlockTridiagonal" begin
        A = BlockTridiagonal(Vcat([fill(1.0,2,1),Matrix(1.0I,2,2),Matrix(1.0I,2,2),Matrix(1.0I,2,2)],Fill(Matrix(1.0I,2,2), ∞)),
                            Vcat([zeros(1,1)], Fill(zeros(2,2), ∞)),
                            Vcat([fill(1.0,1,2),Matrix(1.0I,2,2)], Fill(Matrix(1.0I,2,2), ∞)))

        @test A isa InfiniteLinearAlgebra.BlockTriPertToeplitz
        @test isblockbanded(A)

        @test A[Block.(1:2),Block(1)] == A[1:3,1:1] == reshape([0.,1.,1.],3,1)

        @test BlockBandedMatrix(A)[1:100,1:100] == BlockBandedMatrix(A,(2,1))[1:100,1:100] == BlockBandedMatrix(A,(1,1))[1:100,1:100] == A[1:100,1:100]

        @test (A - I)[1:100,1:100] == A[1:100,1:100]-I
        @test (A + I)[1:100,1:100] == A[1:100,1:100]+I
        @test (I + A)[1:100,1:100] == I+A[1:100,1:100]
        @test (I - A)[1:100,1:100] == I-A[1:100,1:100]

        @test (A - im*I)[1:100,1:100] == A[1:100,1:100]-im*I
        @test (A + im*I)[1:100,1:100] == A[1:100,1:100]+im*I
        @test (im*I + A)[1:100,1:100] == im*I+A[1:100,1:100]
        @test (im*I - A)[1:100,1:100] == im*I-A[1:100,1:100]

        T = mortar(LazyBandedMatrices.Tridiagonal(Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞)))
        #TODO: copy BlockBidiagonal code from BlockBandedMatrices to LazyBandedMatrices
        @test T[Block(2,2)] == [1 2; 3 4]
        @test_broken T[Block(1,3)] == Zeros(2,2)
    end

    @testset "BlockBidiagonal" begin
        B = mortar(LazyBandedMatrices.Bidiagonal(Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞), :U))
        #TODO: copy BlockBidiagonal code from BlockBandedMatrices to LazyBandedMatrices
        @test B[Block(2,3)] == [1 2; 3 4]
        @test_broken B[Block(1,3)] == Zeros(2,2)
    end

    @testset "Fill" begin
        A = _BandedMatrix(Ones(1,∞),ℵ₀,-1,1)
        @test 1.0 .* A isa BandedMatrix{Float64,<:Fill}
        @test Zeros(∞) .* A ≡ Zeros(∞,∞) .* A ≡ A .* Zeros(1,∞) ≡ A .* Zeros(∞,∞) ≡ Zeros(∞,∞)
        @test Ones(∞) .* A isa BandedMatrix{Float64,<:Ones}
        @test A .* Ones(1,∞) isa BandedMatrix{Float64,<:Ones}
        @test 2.0 .* A isa BandedMatrix{Float64,<:Fill}
        @test A .* 2.0 isa BandedMatrix{Float64,<:Fill}
        @test Eye(∞)*A isa BandedMatrix{Float64,<:Ones}
        @test A*Eye(∞) isa BandedMatrix{Float64,<:Ones}

        @test A*A isa BandedMatrix
        @test (A*A)[1:10,1:10] == BandedMatrix(2 => Ones(8))

        Ã = _BandedMatrix(Fill(1,1,∞), ℵ₀, -1,1)
        @test A*Ã isa BandedMatrix
        @test Ã*A isa BandedMatrix
        @test Ã*Ã isa BandedMatrix

        B = _BandedMatrix(Ones(1,10),ℵ₀,-1,1)
        C = _BandedMatrix(Ones(1,10),10,-1,1)
        D = _BandedMatrix(Ones(1,∞),10,-1,1)

        @test (A*B)[1:10,1:10] == (B*C)[1:10,1:10] == (D*A)[1:10,1:10] == D*B == (C*D)[1:10,1:10] == BandedMatrix(2 => Ones(8))
    end

    @testset "Banded Broadcast" begin
        A = _BandedMatrix((1:∞)',ℵ₀,-1,1)
        @test 2.0 .* A isa BandedMatrix{Float64,<:Adjoint}
        @test A .* 2.0 isa BandedMatrix{Float64,<:Adjoint}
        @test Eye(∞)*A isa BandedMatrix{Float64,<:Adjoint}
        @test A*Eye(∞) isa BandedMatrix{Float64,<:Adjoint}
        A = _BandedMatrix(Vcat((1:∞)',Ones(1,∞)),ℵ₀,0,1)
        @test 2.0 .* A isa BandedMatrix
        @test A .* 2.0 isa BandedMatrix
        @test Eye(∞) * A isa BandedMatrix
        @test A * Eye(∞) isa BandedMatrix
        b = 1:∞
        @test BroadcastStyle(typeof(b)) isa LazyArrayStyle{1}
        @test BroadcastStyle(typeof(A)) isa BandedStyle
        @test BroadcastStyle(LazyArrayStyle{1}(), BandedStyle()) isa LazyArrayStyle{2}
        @test BroadcastStyle(LazyArrayStyle{2}(), BandedStyle()) isa LazyArrayStyle{2}
        @test bandwidths(b .* A) == (0,1)

        @test colsupport(b.*A, 1) == 1:1
        @test Base.replace_in_print_matrix(b.*A, 2,1,"0.0") == " ⋅ "
        @test bandwidths(A .* b) == (0,1)
        @test A .* b' isa BroadcastArray
        @test bandwidths(A .* b') == bandwidths(A .* b')
        @test colsupport(A .* b', 3) == 2:3

        A = _BandedMatrix(Ones{Int}(1,∞),ℵ₀,0,0)'
        B = _BandedMatrix((-2:-2:-∞)',ℵ₀ ,-1,1)
        C = Diagonal( 2 ./ (1:2:∞))
        @test bandwidths(A*(B*C)) == (-1,1)
        @test bandwidths((A*B)*C) == (-1,1)

        A = _BandedMatrix(Ones{Int}(1,∞),ℵ₀,0,0)'
        B = _BandedMatrix((-2:-2:-∞)', ℵ₀,-1,1)
        @test MemoryLayout(A+B) isa BroadcastBandedLayout{typeof(+)}
        @test MemoryLayout(2*(A+B)) isa BroadcastBandedLayout{typeof(*)}
        @test bandwidths(A+B) == (0,1)
        @test bandwidths(2*(A+B)) == (0,1)
    end

    @testset "Triangle OP recurrences" begin
        k = mortar(Base.OneTo.(1:∞))
        n = mortar(Fill.(1:∞, 1:∞))
        @test k[Block.(2:3)] isa BlockArray
        @test n[Block.(2:3)] isa BlockArray
        @test k[Block.(2:3)] == [1,2,1,2,3]
        @test n[Block.(2:3)] == [2,2,3,3,3]
        @test blocksize(BroadcastVector(exp,k)) == (ℵ₀,)
        @test BroadcastVector(exp,k)[Block.(2:3)] == exp.([1,2,1,2,3])
        # BroadcastVector(+,k,n)
    end
    # Multivariate OPs Corollary (3)
    # n = 5
    # BlockTridiagonal(Zeros.(1:∞,2:∞),
    #         (n -> Diagonal(((n+2).+(0:n)))/ (2n + 2)).(0:∞),
    #         Zeros.(2:∞,1:∞))

    @testset "KronTrav" begin
        Δ = BandedMatrix(1 => Ones(∞)/2, -1 => Ones(∞))
        A = KronTrav(Δ, Eye(∞))
        @test A[Block(100,101)] isa BandedMatrix
        @test A[Block(100,100)] isa BandedMatrix
        @test A[Block.(1:5), Block.(1:5)] isa BandedBlockBandedMatrix
        B = KronTrav(Eye(∞), Δ)
        @test B[Block(100,101)] isa BandedMatrix
        @test B[Block(100,100)] isa BandedMatrix
        V = view(A+B, Block.(1:5), Block.(1:5))
        @test MemoryLayout(typeof(V)) isa BroadcastBandedBlockBandedLayout{typeof(+)}
        @test arguments(V) == (A[Block.(1:5),Block.(1:5)],B[Block.(1:5),Block.(1:5)])
        @test (A+B)[Block.(1:5), Block.(1:5)] == A[Block.(1:5), Block.(1:5)] + B[Block.(1:5), Block.(1:5)]

        @test blockbandwidths(A+B) == (1,1)
        @test blockbandwidths(2A) == (1,1)
        @test blockbandwidths(2*(A+B)) == (1,1)

        @test subblockbandwidths(A+B) == (1,1)
        @test subblockbandwidths(2A) == (1,1)
        @test subblockbandwidths(2*(A+B)) == (1,1)
    end
end

include("test_hessenbergq.jl")
include("test_infql.jl")
include("test_infqr.jl")
include("test_inful.jl")
include("test_infcholesky.jl")


text:
#############################################################################
#############################################################################
#
# This file implements polynomial GCD
#
#############################################################################
#############################################################################

"""
The extended euclid algorithm for polynomials modulo prime.
"""
function extended_euclid_alg(a::Polynomial, b::Polynomial, prime::Int)
    old_r, r = mod(a,prime), mod(b,prime)
    old_s, s = one(Polynomial), zero(Polynomial)
    old_t, t = zero(Polynomial), one(Polynomial)

    while !iszero(mod(r,prime))
        q = divide(old_r, r)(prime) |> first
        old_r, r = r, mod(old_r - q*r, prime)
        old_s, s = s, mod(old_s - q*s, prime)
        old_t, t = t, mod(old_t - q*t, prime)
    end
    g, s, t = old_r, old_s, old_t
    @assert mod(s*a + t*b - g, prime) == 0
    return g, s, t
end

"""
The GCD of two polynomials modulo prime.
"""
gcd(a::Polynomial, b::Polynomial, prime::Int) = extended_euclid_alg(a,b,prime) |> first


text:
abstract type SimpleNode{T,O} <: AbstractNode{T,O} end

struct SimpleTree{T,O} <: AbstractTree{T,O} end
Base.getindex(::SimpleTree, node::N) where {N<:SimpleNode} = node

nodetype(::SimpleTree{T,O}) where {T,O} = SimpleNode{T,O}

Base.intersect(ray::Ray{T}, n::SimpleNode{T}, tree::SimpleTree{T}) where T =
    intersect(ray, n)

Base.intersect!(intersection::Intersection{T,O}, n::SimpleNode{T,O}, tree::SimpleTree{T,O}) where {T,O} =
    intersect!(intersection, n)

struct BinaryNode{T,O} <: SimpleNode{T,O}
    left::SimpleNode{T,O}
    right::SimpleNode{T,O}

    bbox::AABB{T}
end

function Base.intersect(ray::Ray{T}, n::BinaryNode{T}) where T
    left = intersect(ray, n.left.bbox)
    right = intersect(ray, n.right.bbox)
    (left || right) && (left && intersect(ray, n.left) ||
                      right && intersect(ray, n.right))
end

function Base.intersect!(intersection::Intersection{T,O}, n::BinaryNode{T,O}) where {T,O}
    intersect(intersection.ray, n.bbox) || return
    intersect!(intersection, n.left)
    intersect!(intersection, n.right)
end

add_binary_node!(::SimpleTree{T,O}, left::SimpleNode{T,O}, right::SimpleNode{T,O}) where {T,O} =
    BinaryNode(left, right, AABB(left.bbox, right.bbox))

struct LeafNode{T,O} <: SimpleNode{T,O}
    objs::Vector{O}

    bbox::AABB{T}
end

Base.intersect(ray::Ray{T}, n::LeafNode{T,O}) where {T,O} =
    any(o -> intersect(ray, o), n.objs)

Base.intersect!(intersection::Intersection{T,O}, n::LeafNode{T,O}) where {T,O} =
    foreach(o -> intersect!(intersection, o), n.objs)

add_leaf_node!(::SimpleTree{T,O}, objs::VI, bbox::AABB{T}) where {T,O<:Intersectable{T},VI<:AbstractVector{O}} =
    LeafNode{T,O}(objs[:], bbox)

is_inner(::BinaryNode) = true
is_leaf(::BinaryNode) = false
is_inner(::LeafNode) = false
is_leaf(::LeafNode) = true

objects(n::LeafNode) = n.objs

export SimpleTree


text:
using Test
using SparseIR

@testset "composite.jl" begin
    @testset "Augmented basis with stat = $stat" for stat in (fermion, boson)
        wmax = 2.0
        β = 1e3
        basis = FiniteTempBasis(stat, β, wmax, 1e-6)
        basis_legg = LegendreBasis(stat, β, 2)
        basis_comp = CompositeBasis([basis_legg, basis])

        # G(τ) = c - e^{-τ*pole}/(1 - e^{-β*pole})
        pole = 1.0
        c = 1e-2
        τ_smpl = TauSampling(basis_comp)
        gτ = c .- exp.(-τ_smpl.sampling_points * pole) / (1 - exp(-β * pole))
        gl_from_τ = fit(τ_smpl, gτ)

        gτ_reconst = evaluate(τ_smpl, gl_from_τ)
        @test isapprox(gτ, gτ_reconst; atol=1e-14 * maximum(abs, gτ), rtol=0)
    end
end


text:
using Makie, VideoIO, Dates, Observables

include("testvideo.jl")

# create the test video if it doesn't exist
testvideo = joinpath(tempdir(), "test.mp4")
if !isfile(testvideo)
    createtestvideo(testvideo)
end

# open the video
f = VideoIO.openvideo(testvideo)

# read one image to get its type <-- this might be unecessary since we might know that in advance for all videos
img = Observable(read(f))

# get an estimate for the number of frames, I'm over-shooting here to be sure we have enough space for them
nframes = 2round(Int, VideoIO.get_duration(f.avin.io)/Microsecond(Second(1))*f.framerate)

# prepare the holders for the frames and the time-stamps
frames = Vector{typeof(img[])}() # Vector{tyepof(img[])}(undef, nframes)
timestamps = String[] # Vector{String}(undef, nframes)
sizehint!(frames, nframes)
sizehint!(timestamps, nframes)

# reset the pointer to the begining of the file <-- there might be a more elegant way to achieve this
f = VideoIO.openvideo(testvideo)

# load up the WHOLE video to memory and get the corresponding time-stamps
while !eof(f)
    push!(frames, read(f))
    # note that I'm saving the time stamps as strings, cause that's all we care about right now, might change later
    t = Time(0) + Millisecond(round(Int, 1000gettime(f)))
    push!(timestamps, string(t))
end

# get the actual number of frames (needed actually later)
nframes = length(timestamps)
sizehint!(frames, nframes)

# all the stuff needed for display
pixelaspectratio = VideoIO.aspect_ratio(f)
h = f.height
w = round(typeof(h), f.width * pixelaspectratio)

AbstractPlotting.set_theme!(Attributes(backgroundcolor = :black))

scene = Scene(resolution = (w, h))
Makie.image!(scene, 1:h, 1:w, img, show_axis = false, scale_plot = false)
Makie.rotate!(scene, -0.5π)
update_limits!(scene)
update_cam!(scene)
update!(scene)

# everything starts at this observable, the current frame
currentframe = Observable(1)

# it controls which frame is displayed
lift(currentframe) do i
    img[] = frames[i]
end

# the slider, updating currentframe and getting its string from timestamps
slider_h = slider!(rsc(), 1:length(timestamps), start = 1, valueprinter = i -> timestamps[i], textcolor = :white)
lift(slider_h[end][:value]) do i
    currentframe[] = i
end

# I have to have a forward button becuase apparently the stuff in the lift gets evaluated, so to have the movie start in frame #1 I need to go forward once and the backwards...
fwdbutton = button!(rsc(), ">", textcolor = :white)
lift(fwdbutton[end][:clicks]) do _
    if currentframe[] < nframes # check we are not in the last frame
        currentframe[] += 1
        slider_h[end][:value][] += 1
        # so now the movie is in frame #2
    end
    0
end
bckbutton = button!(rsc(), "<", textcolor = :white)
lift(bckbutton[end][:clicks]) do _
    if currentframe[] > 1 # check we're not in frame #1
        currentframe[] -= 1
        slider_h[end][:value][] -= 1
        # and now we're back in frame #1
    end
    0
end

# done!
hbox(vbox(bckbutton, slider_h, fwdbutton), scene)


text:
module HSARuntime

# copied from CUDAdrv/src/CUDAdrv.jl
const ext = joinpath(dirname(@__DIR__), "deps", "ext.jl")
isfile(ext) || error("HSARuntime.jl has not been built, please run Pkg.build(\"HSARuntime\").")
include(ext)
if !configured
    # default (non-functional) values for critical variables,
    # making it possible to _load_ the package at all times.
    const libhsaruntime_version = v"0.0"
    const libhsaruntime_vendor = "none"
    const libhsaruntime_path = nothing
end

atexit() do
    configured && hsa_shut_down()
end

### Imports ###

using CEnum
using Setfield
using Libdl

### Exports ###

export HSAAgent, HSAQueue, HSAExecutable, HSAKernelInstance, HSAArray, HSASignal
export get_agents, profile, get_first_isa, launch!
export get_default_agent, get_default_queue

### HSA Runtime Wrapper ###

include(joinpath(@__DIR__, "..", "deps", "hsa_types.jl"))
include(joinpath(@__DIR__, "..", "deps", "hsa_interface.jl"))
include(joinpath(@__DIR__, "..", "deps", "hsa_extras.jl"))

### HSA Errors ###

export HSAError

struct HSAError <: Exception
    code::hsa_status_t
end

"""
    description(err::HSAError)

Gets the string description of an error code.
"""
function description(err::HSAError)
    str_ref = Ref{Cstring}()
    hsa_status_string(err.code, str_ref)
    unsafe_string(str_ref[])
end

function Base.showerror(io::IO, err::HSAError)
    print(io, "HSA error (code #$(Int(err.code)), $(description(err)))")
end

Base.show(io::IO, err::HSAError) = print(io, "HSAError($(Int(err.code)), $(description(err)))")

function check(result)
    if result != HSA_STATUS_SUCCESS
        throw(HSAError(result))
    end
end

### Types ###

mutable struct HSAAgent
    agent::hsa_agent_t
end

mutable struct HSAQueue
    agent::HSAAgent
    queue::Ref{Ptr{hsa_queue_t}}
end

mutable struct HSAExecutable
    agent::HSAAgent
    executable::Ref{hsa_executable_t}
    code_object::Ref{hsa_code_object_t}
    data::Vector{UInt8}
end

mutable struct HSAKernelInstance{T<:Tuple}
    agent::HSAAgent
    exe::HSAExecutable
    sym::String
    args::T
    # FIXME: Stop being lazy and put in some damn types!
    kernel_object
    kernarg_segment_size
    group_segment_size
    private_segment_size
    kernarg_address::Ref{Ptr{Nothing}}
end

mutable struct HSAArray{T,N} <: AbstractArray{T,N}
    size::Dims{N}
    handle::Ptr{T}
end

mutable struct HSASignal
    signal::Ref{hsa_signal_t}
end

### Constants ###

const DEFAULT_AGENT = Ref{HSAAgent}()
const DEFAULT_QUEUES = IdDict{HSAAgent,HSAQueue}()

### @cfunction Callbacks ###

function agent_iterate_isas_cb(isa::hsa_isa_t, ptr::Ptr{hsa_isa_t})
    unsafe_store!(ptr, isa)
    return HSA_STATUS_INFO_BREAK
end

function iterate_agents_cb(agent::hsa_agent_t, agents)
    push!(agents, HSAAgent(agent))
    return HSA_STATUS_SUCCESS
end

function finalizer_iterate_isas_cb(isa::hsa_isa_t, data)
    @info "Finalizer ISAs:"

    name = repeat(" ", 64)
    hsa_isa_get_info_alt(isa, HSA_ISA_INFO_NAME, name) |> check
    @info "    ISA Name: $name"

    return HSA_STATUS_SUCCESS
end

"Determines if a memory region can be used for kernarg allocations."
function get_kernarg_memory_region_cb(region::hsa_region_t, data::Ptr{hsa_region_t})
    segment = Ref{hsa_region_segment_t}()
    hsa_region_get_info(region, HSA_REGION_INFO_SEGMENT, segment)
    if (segment[] != HSA_REGION_SEGMENT_GLOBAL)
        return HSA_STATUS_SUCCESS
    end

    flags = Ref{hsa_region_global_flag_t}()
    hsa_region_get_info(region, HSA_REGION_INFO_GLOBAL_FLAGS, flags)
    if (flags[] & HSA_REGION_GLOBAL_FLAG_KERNARG > 0)
        unsafe_store!(data, region)
        return HSA_STATUS_INFO_BREAK
    end

    return HSA_STATUS_SUCCESS
end

"Determines if a memory region can be used for fine grained allocations."
function get_fine_grained_memory_region_cb(region::hsa_region_t, data::Ptr{hsa_region_t})
    segment = Ref{hsa_region_segment_t}()
    hsa_region_get_info(region, HSA_REGION_INFO_SEGMENT, segment)
    if (segment[] != HSA_REGION_SEGMENT_GLOBAL)
        return HSA_STATUS_SUCCESS
    end

    flags = Ref{hsa_region_global_flag_t}()
    hsa_region_get_info(region, HSA_REGION_INFO_GLOBAL_FLAGS, flags)
    if (flags[] & HSA_REGION_GLOBAL_FLAG_FINE_GRAINED > 0)
        unsafe_store!(data, region)
        return HSA_STATUS_INFO_BREAK
    end

    return HSA_STATUS_SUCCESS
end

function iterate_exec_prog_syms_cb(exe::hsa_executable_t, agent::hsa_agent_t, sym::hsa_executable_symbol_t, data)
    name = repeat(" ", 64)
    hsa_executable_symbol_get_info(sym, HSA_EXECUTABLE_SYMBOL_INFO_NAME, name) |> check
    @info "    Symbol Name: $name"

    return HSA_STATUS_SUCCESS
end

### Constructors and Finalizers ###

function HSAQueue(agent::HSAAgent)
    queue_size = Ref{UInt32}(0)
    hsa_agent_get_info(agent.agent, HSA_AGENT_INFO_QUEUE_MAX_SIZE, queue_size) |> check
    @assert queue_size[] > 0
    queue = HSAQueue(agent, Ref{Ptr{hsa_queue_t}}())
    hsa_queue_create(agent.agent, queue_size[], HSA_QUEUE_TYPE_SINGLE,
            C_NULL, C_NULL, typemax(UInt32), typemax(UInt32), queue.queue) |> check
    finalizer(queue) do queue
        hsa_queue_destroy(queue.queue[]) |> check
    end
    return queue
end
function HSASignal()
    signal = HSASignal(Ref{hsa_signal_t}())
    hsa_signal_create(1, 0, C_NULL, signal.signal)
    finalizer(signal) do signal
        hsa_signal_destroy(signal.signal[]) |> check
    end
    return signal
end
function HSAExecutable(agent::HSAAgent, data::Vector{UInt8}, symbol::String)
    #= NOTE: Everything I can see indicates that profile is always FULL
    profile = Ref{hsa_profile_t}()
    hsa_agent_get_info(agent.agent, HSA_AGENT_INFO_PROFILE, profile) |> check
    =#

    code_object = Ref{hsa_code_object_t}(hsa_code_object_t(0))
    hsa_code_object_deserialize(data, sizeof(data), C_NULL, code_object) |> check
    @assert code_object[].handle != 0

    executable = Ref{hsa_executable_t}()
    hsa_executable_create(HSA_PROFILE_FULL, HSA_EXECUTABLE_STATE_UNFROZEN, C_NULL, executable) |> check
    hsa_executable_load_code_object(executable[], agent.agent, code_object[], C_NULL) |> check
    hsa_executable_freeze(executable[], "") |> check

    exe = HSAExecutable(agent, executable, code_object, data)
    finalizer(exe) do exe
        hsa_executable_destroy(exe.executable[]) |> check
        hsa_code_object_destroy(exe.code_object[]) |> check
    end
    return exe
end

function HSAKernelInstance(agent::HSAAgent, exe::HSAExecutable, symbol::String, args::Tuple)
    #= TODO: Make this available for debugging purposes
    func = @cfunction(iterate_exec_prog_syms_cb, hsa_status_t,
            (hsa_executable_t, hsa_agent_t, hsa_executable_symbol_t, Ptr{Cvoid}))
    hsa_executable_iterate_agent_symbols(exe.executable[], agent.agent, func, C_NULL) |> check
    =#

    exec_symbol = Ref{hsa_executable_symbol_t}()
    hsa_executable_get_symbol(exe.executable[], C_NULL, symbol, agent.agent, 0, exec_symbol) |> check

    kernel_object = Ref{UInt64}(0)
    hsa_executable_symbol_get_info(exec_symbol[], HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_OBJECT, kernel_object) |> check

    kernarg_segment_size = Ref{UInt32}(0)
    hsa_executable_symbol_get_info(exec_symbol[], HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_KERNARG_SEGMENT_SIZE, kernarg_segment_size) |> check

    group_segment_size = Ref{UInt32}(0)
    hsa_executable_symbol_get_info(exec_symbol[], HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_GROUP_SEGMENT_SIZE, group_segment_size) |> check

    private_segment_size = Ref{UInt32}(0)
    hsa_executable_symbol_get_info(exec_symbol[], HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_PRIVATE_SEGMENT_SIZE, private_segment_size) |> check

    finegrained_region = get_region(agent, :finegrained)
    kernarg_region = get_region(agent, :kernarg)

    # Allocate the kernel argument buffer from the correct region
    kernarg_address = Ref{Ptr{Nothing}}()
    hsa_memory_allocate(kernarg_region[], kernarg_segment_size[], kernarg_address) |> check
    ctr = 0x0
    for arg in args
        rarg = Ref(arg)
        ccall(:memcpy, Cvoid,
            (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
            kernarg_address[]+ctr, rarg, sizeof(arg))
        ctr += sizeof(arg)
    end

    kernel = HSAKernelInstance(agent, exe, symbol, args, kernel_object,
                               kernarg_segment_size, group_segment_size,
                               private_segment_size, kernarg_address)
    finalizer(kernel) do kernel
        hsa_memory_free(kernel.kernarg_address[]) |> check
    end
    return kernel
end

# TODO: Support non-isbitstype allocations
function HSAArray(agent::HSAAgent, ::Type{T}, size::NTuple{N,Int}) where {T,N}
    @assert isprimitivetype(T) "$T is not a primitive type"
    @assert all(x->x>0, size) "Invalid array size: $size"
    region = get_region(agent, :finegrained)
    nbytes = sizeof(T) * prod(size)
    handle = Ref{Ptr{T}}()
    hsa_memory_allocate(region[], nbytes, handle) |> check
    arr = HSAArray{T,N}(size, handle[])
    finalizer(arr) do arr
        hsa_memory_free(arr.handle) |> check
    end
    return arr
end
HSAArray(::Type{T}, size::NTuple{N,Int}) where {T,N} =
    HSAArray(DEFAULT_AGENT[], T, size)

### HSAArray Methods ###

function HSAArray(agent::HSAAgent, arr::Array{T,N}) where {T,N}
    harr = HSAArray(agent, T, size(arr))
    for idx in eachindex(arr)
        harr[idx] = arr[idx]
    end
    return harr
end
HSAArray(arr::Array{T,N}) where {T,N} =
    HSAArray(DEFAULT_AGENT[], arr)

function Array(harr::HSAArray{T,N}) where {T,N}
    arr = Array{T}(undef, size(harr))
    ref_arr = Ref(arr)
    GC.@preserve ref_arr begin
        ccall(:memcpy, Cvoid,
            (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
            ref_arr, harr.handle, sizeof(arr))
    end
    return harr
end

Base.pointer(arr::HSAArray) = arr.handle
Base.IndexStyle(::Type{<:HSAArray}) = Base.IndexLinear()
Base.IndexStyle(::HSAArray) = Base.IndexLinear()
function Base.iterate(A::HSAArray, i=1) # copy-pasta from Base
    Base.@_inline_meta
    (i % UInt) - 1 < length(A) ? (@inbounds A[i], i + 1) : nothing
end
Base.similar(arr::HSAArray{T,N}) where {T,N} =
    HSAArray(T, size(arr))
Base.similar(agent, arr::HSAArray{T,N}) where {T,N} =
    HSAArray(agent, T, size(arr))

Base.size(arr::HSAArray) = arr.size
Base.length(arr::HSAArray) = sum(size(arr))

function Base.fill!(arr::HSAArray{T,N}, value::T) where {T,N}
    for idx in 1:length(arr)
        arr[idx] = value
    end
end
@inline function Base.getindex(arr::HSAArray{T,N}, idx) where {T,N}
    @boundscheck checkbounds(arr, idx)
    Base.unsafe_load(pointer(arr), idx)::T
end
@inline function Base.setindex!(arr::HSAArray{T,N}, value, idx) where {T,N}
    @boundscheck checkbounds(arr, idx)
    Base.unsafe_store!(pointer(arr), value, idx)
end

### Methods ###

function Base.show(io::IO, agent::HSAAgent)
    print(io, "HSAAgent($(agent.agent)): Name=$(get_name(agent)), Type=$(device_type(agent))")
end

function get_agents()
    agents = Ref(Vector{HSAAgent}())
    GC.@preserve agents begin
        func = @cfunction(iterate_agents_cb, hsa_status_t,
                (hsa_agent_t, Ref{Vector{HSAAgent}}))
        hsa_iterate_agents(func, agents) |> check
        _agents = agents[]
    end
    return _agents
end
get_agents(kind::Symbol) =
    filter(agent->device_type(agent)==kind, get_agents())

get_default_agent() = DEFAULT_AGENT[]
function set_default_agent!(kind::Symbol)
    DEFAULT_AGENT[] = first(get_agents(kind))
end
set_default_agent!() = set_default_agent!(:gpu)

"Gets the default queue for an agent, creating it if necessary"
function get_default_queue(agent::HSAAgent)
    queue = get!(DEFAULT_QUEUES, agent, HSAQueue(agent))
    return queue
end
get_default_queue() =
    get_default_queue(get_default_agent())

function get_name(agent::HSAAgent)
    # TODO: Get name length first!
    name = repeat(" ", 64)
    hsa_agent_get_info(agent.agent, HSA_AGENT_INFO_NAME, name) |> check
    return name
end
function profile(agent::HSAAgent)
    profile = Ref{hsa_profile_t}()
    hsa_agent_get_info(agent.agent, HSA_AGENT_INFO_PROFILE, profile) |> check
    return profile[]
end

function device_type(agent::HSAAgent)
    devtype = Ref{UInt32}()
    hsa_agent_get_info(agent.agent, HSA_AGENT_INFO_DEVICE, devtype) |> check
    if hsa_device_type_t(devtype[]) == HSA_DEVICE_TYPE_CPU
        return :cpu
    elseif hsa_device_type_t(devtype[]) == HSA_DEVICE_TYPE_GPU
        return :gpu
    elseif hsa_device_type_t(devtype[]) == HSA_DEVICE_TYPE_DSP
        return :dsp
    end
end

# TODO: Get all ISAs
function get_first_isa(agent::HSAAgent)
    isa = Ref{hsa_isa_t}()
    func = @cfunction(agent_iterate_isas_cb, hsa_status_t,
        (hsa_isa_t, Ptr{hsa_isa_t}))
    ret = hsa_agent_iterate_isas(agent.agent, func, isa)
    @assert ret == HSA_STATUS_INFO_BREAK "Failed to find an agent ISA"

    len = Ref{Cuint}(0)
    hsa_isa_get_info_alt(isa[], HSA_ISA_INFO_NAME_LENGTH, len) |> check
    name = repeat(" ", len[])
    hsa_isa_get_info_alt(isa[], HSA_ISA_INFO_NAME, name) |> check
    isa_name = string(rstrip(last(split(name, "-")), '\0'))
    return isa_name
end

function get_region(agent::HSAAgent, kind::Symbol)
    region = newref!(Ref{hsa_region_t}, typemax(UInt64))
    if kind == :finegrained
        func = @cfunction(get_fine_grained_memory_region_cb, hsa_status_t,
                (hsa_region_t, Ptr{hsa_region_t}))
    elseif kind == :kernarg
        func = @cfunction(get_kernarg_memory_region_cb, hsa_status_t,
                (hsa_region_t, Ptr{hsa_region_t}))
    else
        error("Region kind $kind not supported")
    end
    hsa_agent_iterate_regions(agent.agent, func, region)
    check((region[].handle == typemax(UInt64)) ? HSA_STATUS_ERROR : HSA_STATUS_SUCCESS)
    return region
end

function launch!(queue::HSAQueue, kernel::HSAKernelInstance, signal::HSASignal;
                 workgroup_size=nothing, grid_size=nothing)
    @assert workgroup_size !== nothing "Must specify workgroup_size kwarg"
    @assert grid_size !== nothing "Must specify grid_size kwarg"

    queue = queue.queue
    signal = signal.signal
    args = kernel.args

    # Obtain the current queue write index
    index = hsa_queue_load_write_index_relaxed(queue[])

    # TODO: Make this less ugly
    dispatch_packet = Ref{hsa_kernel_dispatch_packet_t}()
    ccall(:memset, Cvoid,
        (Ptr{Cvoid}, Cint, Csize_t),
        dispatch_packet, 0, sizeof(dispatch_packet[]))
    _packet = dispatch_packet[]
    @set! _packet.setup = 0
    @set! _packet.setup |= 1 << Int(HSA_KERNEL_DISPATCH_PACKET_SETUP_DIMENSIONS)
    @set! _packet.workgroup_size_x = workgroup_size[1]
    @set! _packet.workgroup_size_y = workgroup_size[2]
    @set! _packet.workgroup_size_z = workgroup_size[3]
    @set! _packet.grid_size_x = grid_size[1]
    @set! _packet.grid_size_y = grid_size[2]
    @set! _packet.grid_size_z = grid_size[3]
    @set! _packet.completion_signal = signal[]
    @set! _packet.kernel_object = kernel.kernel_object[]
    @set! _packet.kernarg_address = kernel.kernarg_address[]
    @set! _packet.private_segment_size = kernel.private_segment_size[]
    @set! _packet.group_segment_size = kernel.group_segment_size[]
    dispatch_packet = Ref{hsa_kernel_dispatch_packet_t}(_packet)

    _queue = unsafe_load(queue[])
    queueMask = UInt32(_queue.size - 1)
    baseaddr_ptr = Ptr{hsa_kernel_dispatch_packet_t}(_queue.base_address)
    baseaddr_ptr += sizeof(hsa_kernel_dispatch_packet_t) * (index & queueMask)
    dispatch_packet_ptr = Base.unsafe_convert(Ptr{hsa_kernel_dispatch_packet_t}, dispatch_packet)
    unsafe_copyto!(baseaddr_ptr, dispatch_packet_ptr, 1)

    # Atomically store the header
    atomic_store_n!(x::Ptr{UInt16}, v::UInt16) =
        Base.llvmcall("""
        %ptr = inttoptr i$(Sys.WORD_SIZE) %0 to i16*
        store atomic i16 %1, i16* %ptr release, align 16
        ret void
        """, Cvoid, Tuple{Ptr{UInt16}, UInt16},
        Base.unsafe_convert(Ptr{UInt16}, x), v)
    header = Ref{UInt16}(0)
    header[] |= Int(HSA_FENCE_SCOPE_SYSTEM) << Int(HSA_PACKET_HEADER_ACQUIRE_FENCE_SCOPE)
    header[] |= Int(HSA_FENCE_SCOPE_SYSTEM) << Int(HSA_PACKET_HEADER_RELEASE_FENCE_SCOPE)
    header[] |= Int(HSA_PACKET_TYPE_KERNEL_DISPATCH) << Int(HSA_PACKET_HEADER_TYPE)
    atomic_store_n!(Base.unsafe_convert(Ptr{UInt16}, baseaddr_ptr), header[])

    # Increment the write index and ring the doorbell to dispatch the kernel
    hsa_queue_store_write_index_relaxed(queue[], index+1)
    hsa_signal_store_relaxed(_queue.doorbell_signal, index)
end

# FIXME: Make this non-blocking for libuv
function Base.wait(signal::HSASignal)
    # Wait on the dispatch completion signal until the kernel is finished
    value = Ref{hsa_signal_value_t}()
    value[] = hsa_signal_wait_acquire(signal.signal[],
                                      HSA_SIGNAL_CONDITION_LT, 1, typemax(UInt64),
                                      HSA_WAIT_STATE_BLOCKED)
end

include("memory.jl")

function __init__()
    # We want to always be able to load the package, regardless of
    # whether HSA libraries are found (just like the CUDA* packages)
    configured || return

    # Make sure we load the library found by the last `] build`
    push!(Libdl.DL_LOAD_PATH, dirname(libhsaruntime_path))
    # TODO: Do the same (if possible) for the debug library

    hsa_init() |> check

    agents = get_agents(:gpu)
    if length(agents) > 0
        DEFAULT_AGENT[] = first(agents)
    end
end

end


text:
# Generate documentation with this command:
# (cd docs && julia --color=yes make.jl)

push!(LOAD_PATH, "..")

using Documenter
using AbstractSphericalHarmonics

makedocs(; sitename="AbstractSphericalHarmonics", format=Documenter.HTML(), modules=[AbstractSphericalHarmonics])

deploydocs(; repo="github.com/eschnett/AbstractSphericalHarmonics.jl.git", devbranch="main", push_preview=true)


text:
using BinaryBuilder

# Collection of pre-build pandoc binaries
name = "pandoc"
version = v"2.14"
pandoc_ver = "2.14"

url_prefix = "https://github.com/jgm/pandoc/releases/download/$(pandoc_ver)/pandoc-$(pandoc_ver)"
sources = [
    ArchiveSource("$(url_prefix)-linux-amd64.tar.gz", "e24cca7d241d907d690428a296861cbcbdaf4686ef887ac65581a930b2e6038c"; unpack_target = "x86_64-linux-gnu"),
    ArchiveSource("$(url_prefix)-macOS.zip", "630564487c30602ee299b463efc4fb48b418a6f254c123a11f3748426b2007b3"; unpack_target = "x86_64-apple-darwin14"),
    ArchiveSource("$(url_prefix)-windows-x86_64.zip", "341935e8de2bc58153b6aabcced8bf73d753bc6b97efdf0128f52da82841b0c5"; unpack_target = "x86_64-w64-mingw32"),
    ArchiveSource("$(url_prefix)-linux-arm64.tar.gz", "fd334f48c12b17ec264b2c49a554c584cdb29828c29fbbde7c1f57232a30975d"; unpack_target = "aarch64-linux-gnu"),
    FileSource("https://raw.githubusercontent.com/jgm/pandoc/$(pandoc_ver)/COPYRIGHT", "adcfa50add0dd23fda4937830ee6401a45638cffa15b9b33b1932f833f4fab75"),
    FileSource("https://raw.githubusercontent.com/jgm/pandoc/$(pandoc_ver)/COPYING.md", "e7ea3adeab955103a837b692ca0017cb3abbed0d3dccbfa499d6b2b825d698c3"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
mkdir -p "${bindir}"
if [[ "${target}" != *-mingw* ]]; then
    subdir="bin/"
fi
cp ${target}/pandoc-*/${subdir}pandoc${exeext} ${bindir}
chmod +x ${bindir}/*
install_license COPYRIGHT
install_license COPYING.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
    Platform("aarch64", "linux")
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("pandoc", :pandoc),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    PBVISolver <: Solver

Options dictionary for Point-Based Value Iteration for POMDPs.

# Fields
- `max_iterations::Int64` the maximal number of iterations the solver runs. Default: 10
- `ϵ::Float64` the maximal gap between alpha vector improve steps. Default = 0.01
- `verbose::Bool` switch for solver text output. Default: false
"""
struct PBVISolver <: Solver
    max_iterations::Int64
    ϵ::Float64
    verbose::Bool
end

function PBVISolver(;max_iterations::Int64=10, ϵ::Float64=0.01, verbose::Bool=false)
    return PBVISolver(max_iterations, ϵ, verbose)
end

"""
    AlphaVec

Pair of alpha vector and corresponding action.

# Fields
- `alpha` α vector
- `action` action corresponding to α vector
"""
struct AlphaVec
    alpha::Vector{Float64}
    action::Any
end

==(a::AlphaVec, b::AlphaVec) = (a.alpha,a.action) == (b.alpha, b.action)
Base.hash(a::AlphaVec, h::UInt) = hash(a.alpha, hash(a.action, h))

function _argmax(f, X)
    return X[argmax(map(f, X))]
end

# adds probabilities of terminals in b to b′ and normalizes b′
function belief_norm(pomdp, b, b′, terminals, not_terminals)
    if sum(b′[not_terminals]) != 0.
        if !isempty(terminals)
            b′[not_terminals] = b′[not_terminals] / (sum(b′[not_terminals]) / (1. - sum(b[terminals]) - sum(b′[terminals])))
            b′[terminals] += b[terminals]
        else
            b′[not_terminals] /= sum(b′[not_terminals])
        end
    else
        b′[terminals] += b[terminals]
        b′[terminals] /= sum(b′[terminals])
    end
    return b′
end

# Backups belief with α vector maximizing dot product of itself with belief b
function backup_belief(pomdp::POMDP, Γ, b)
    S = ordered_states(pomdp)
    A = ordered_actions(pomdp)
    O = ordered_observations(pomdp)
    γ = discount(pomdp)
    r = StateActionReward(pomdp)

    Γa = Vector{Vector{Float64}}(undef, length(A))

    not_terminals = [stateindex(pomdp, s) for s in S if !isterminal(pomdp, s)]
    terminals = [stateindex(pomdp, s) for s in S if isterminal(pomdp, s)]
    for a in A
        Γao = Vector{Vector{Float64}}(undef, length(O))
        trans_probs = dropdims(sum([pdf(transition(pomdp, S[is], a), sp) * b.b[is] for sp in S, is in not_terminals], dims=2), dims=2)
        if !isempty(terminals) trans_probs[terminals] .+= b.b[terminals] end

        for o in O
            # update beliefs
            obs_probs = pdf.(map(sp -> observation(pomdp, a, sp), S), [o])
            b′ = obs_probs .* trans_probs

            if sum(b′) > 0.
                b′ = DiscreteBelief(pomdp, b.state_list, belief_norm(pomdp, b.b, b′, terminals, not_terminals))
            else
                b′ = DiscreteBelief(pomdp, b.state_list, zeros(length(S)))
            end

            # extract optimal alpha vector at resulting belief
            Γao[obsindex(pomdp, o)] = _argmax(α -> α ⋅ b′.b, Γ)
        end

        # construct new alpha vectors
        Γa[actionindex(pomdp, a)] = [r(s, a) + (!isterminal(pomdp, s) ? (γ * sum(pdf(transition(pomdp, s, a), sp) * pdf(observation(pomdp, s, a, sp), o) * Γao[i][j]
                                        for (j, sp) in enumerate(S), (i, o) in enumerate(O))) : 0.)
                                        for s in S]
    end

    # find the optimal alpha vector
    idx = argmax(map(αa -> αa ⋅ b.b, Γa))
    alphavec = AlphaVec(Γa[idx], A[idx])

    return alphavec
end

# Iteratively improves α vectors until the gap between steps is lesser than ϵ
function improve(pomdp, B, Γ, solver)
    alphavecs = nothing
    while true
        Γold = Γ
        alphavecs = [backup_belief(pomdp, Γold, b) for b in B]
        Γ = [alphavec.alpha for alphavec in alphavecs]
        prec = max([sum(abs.(dot(α1, b.b) .- dot(α2, b.b))) for (α1, α2, b) in zip(Γold, Γ, B)]...)
        if solver.verbose println("    Improving alphas, maximum gap between old and new α vector: $(prec)") end
        prec > solver.ϵ || break
    end

    return Γ, alphavecs
end

# Returns all possible, not yet visited successors of current belief b
function successors(pomdp, b, Bs)
    S = ordered_states(pomdp)
    not_terminals = [stateindex(pomdp, s) for s in S if !isterminal(pomdp, s)]
    terminals = [stateindex(pomdp, s) for s in S if isterminal(pomdp, s)]
    succs = []

    for a in actions(pomdp)
        trans_probs = dropdims(sum([pdf(transition(pomdp, S[is], a), sp) * b[is] for sp in S, is in not_terminals], dims=2), dims=2)
        if !isempty(terminals) trans_probs[terminals] .+= b[terminals] end

        for o in observations(pomdp)
            #update belief
            obs_probs = pdf.(map(sp -> observation(pomdp, a, sp), S), [o])
            b′ = obs_probs .* trans_probs

            if sum(b′) > 0.
                b′ = belief_norm(pomdp, b, b′, terminals, not_terminals)

                if !in(b′, Bs)
                    push!(succs, b′)
                end
            end
        end
    end

    return succs
end

# Computes distance of successor to the belief vectors in belief space
function succ_dist(pomdp, bp, B)
    dist = [norm(bp - b.b, 1) for b in B]
    return max(dist...)
end

# Expands the belief space with the most distant belief vector
# Returns new belief space, set of belifs and early termination flag
function expand(pomdp, B, Bs)
    B_new = copy(B)
    for b in B
        succs = successors(pomdp, b.b, Bs)
        if length(succs) > 0
            b′ = succs[argmax([succ_dist(pomdp, bp, B) for bp in succs])]
            push!(B_new, DiscreteBelief(pomdp, b′))
            push!(Bs, b′)
        end
    end

    return B_new, Bs, length(B) == length(B_new)
end

# 1: B ← {b0}
# 2: while V has not converged to V∗ do
# 3:    Improve(V, B)
# 4:    B ← Expand(B)
function solve(solver::PBVISolver, pomdp::POMDP)
    S = ordered_states(pomdp)
    A = ordered_actions(pomdp)
    γ = discount(pomdp)
    r = StateActionReward(pomdp)

    # best action worst state lower bound
    α_init = 1 / (1 - γ) * maximum(minimum(r(s, a) for s in S) for a in A)
    Γ = [fill(α_init, length(S)) for a in A]

    #init belief, if given distribution, convert to vector
    init = initialize_belief(DiscreteUpdater(pomdp), initialstate(pomdp))
    B = [init]
    Bs = Set([init.b])

    if solver.verbose println("Running PBVI solver on $(typeof(pomdp)) problem with following settings:\n    max_iterations = $(solver.max_iterations), ϵ = $(solver.ϵ), verbose = $(solver.verbose)\n+----------------------------------------------------------+") end

    # original code should run until V converges to V*, this yet needs to be implemented
    # for example as: while max(@. abs(newV - oldV)...) > solver.ϵ
    # However this probably would not work, as newV and oldV have different number of elements (arrays of alphas)
    alphavecs = nothing
    for i in 1:solver.max_iterations
        Γ, alphavecs = improve(pomdp, B, Γ, solver)
        B, Bs, early_term = expand(pomdp, B, Bs)
        if solver.verbose println("Iteration $(i) executed, belief set contains $(length(Bs)) belief vectors.") end
        if early_term
            if solver.verbose println("Belief space did not expand. \nTerminating early.") end
            break
        end
    end

    if solver.verbose println("+----------------------------------------------------------+") end
    acts = [alphavec.action for alphavec in alphavecs]
    return AlphaVectorPolicy(pomdp, Γ, acts)
end

@POMDPLinter.POMDP_require solve(solver::PBVISolver, pomdp::POMDP) begin
    P = typeof(pomdp)
    S = state_type(P)
    A = action_type(P)
    O = observation_type(P)
    @req discount(::P) # discount factor
    @subreq ordered_states(pomdp)
    @subreq ordered_actions(pomdp)
    @subreq ordered_observations(pomdp)
    @req transition(::P,::S,::A)
    @req reward(::P,::S,::A)
    ss = states(pomdp)
    as = actions(pomdp)
    os = observations(pomdp)
    @req length(::typeof(ss))
    s = first(iterator(ss))
    a = first(iterator(as))
    dist = transition(pomdp, s, a)
    D = typeof(dist)
    @req pdf(::D,::S)

    odist = observation(pomdp, a, s)
    OD = typeof(odist)
    @req pdf(::OD,::O)
end


text:
using ArgParse
using BenchmarkTools

include("configurations\\MasterConfiguration.jl")
include("master\\ODP_MasterAlgorithm.jl")
include("utils\\Result.jl")

include("utils\\PlotsUtils.jl")

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "n"
            help = "n dededede"
            arg_type = Int
            required = true
        "k"
            help = "k dededede"
            arg_type = Int
            required = true
        "m"
            help = "m dededede"
            arg_type = Int
            required = true
        "port"
            help = "Server port"
            arg_type = Int
            default = 10000
        "dist"
            help = "Distance function"
            arg_type = String
            default = "square"
        "conn_type"
            help = "Connection type"
            arg_type = String
            default = "tcp"
    end
    return parse_args(s)
end
function main()
    # @show parsed_args = parse_commandline()
    # println("Parsed args:")
    # for (arg,val) in parsed_args
    #     print("  $arg  =>  ")
    #     show(val)
    #     println()
    # end

    #  Setting up Slave
    filepath1 = "datasets\\make_moons100k.txt"
    filepath2 = "datasets\\make_moons100k.txt"
    filepath3 = "datasets\\make_moons100k.txt"
    filepath4 = "datasets\\make_moons100k.txt"

    generateNewFile = false

    args_conf_slave = [filepath1, generateNewFile]

    if args_conf_slave[2]
        args_conf_slave[1] = "datasets\\test3.txt"
        rndMatrix = map(x->300000.0+x*(1000),rand(10000,2))
        rndMatrix3 = map(x->3000.0+x*(1005),rand(100,2))
        rndMatrix2 = map(x->3000.0+x*(1002),rand(100,2))

        foo = [rndMatrix, rndMatrix2, rndMatrix3]
        A = reduce(vcat,foo)

        generatePoints(filepath1, A)
    end

    filepath1 = args_conf_slave[1]

    #  Setting up Master
    args_conf_master = [20,6,6,"euclidean"]
    n = args_conf_master[1]
    k = args_conf_master[2]
    m = args_conf_master[3]
    dist = args_conf_master[4]

    config1 = MasterConfiguration(100, 20, 20, dist)

    masterAlg1 = getMasterAlgorithm(config1)

    result1 = runODPAlgorithm(masterAlg1, filepath1)
    result2 = runODPAlgorithm(masterAlg1, filepath2)
    result3 = runODPAlgorithm(masterAlg1, filepath3)
    result4 = runODPAlgorithm(masterAlg1, filepath4)

    r = Result[result1,result2,result3,result4]

    return r

end

r = main()

# Print report..
print(getReport(r))

# Plotting..
display(plotSets("Dataset", r[4].dataset, r[4].solution, ["Dataset","solution"]))

c1 = 1
for i in r
    display(plotSets("Dataset"*string(count), i.dataset, i.solvingSet, i.solution, ["Dataset","SolvingSet","solution"]))
    global c1+=1
end

#-----------------------------------------------------------
# ()|D|, |S|) - plot
x = Vector{Int64}()
y1 = Vector{Int64}()
for i in r
    push!(x,i.dataset_size)
    push!(y1,i.solvingSet_size)
end
display(plot(x, y1, title = "Scaling analysis of the ODP solving sets sizes", label = "|S| -> n=100 k=20", marker=[:o], c=[:orange], legend=:right))

# (|D|, execution time in seconds) - plot
y2 = [21.04, 39.5, 90.45,207.5]
display(plot(x, y2, title = "Scaling analysis of the execution times", label = "time(s) -> n=100 k=20", marker=[:o], c=[:red], legend=:right))

# for i in 1:length(result1.itC)
#     plt = plotSets("Iteration#"*string(i), result1.dataset, result1.itS[i], result1.itC[i], ["Dataset","SolvingSet","Cand"])
#     display(plt)
# end
# plotSets(result1.dataset, result1.solvingSet, result1.solution, ["Dataset","SolvingSet","solution"])
#
# plotSets("Full Dataset",result1.dataset, ["Dataset"])
#plotSets("Solution", result1.solvingSet, result1.solution, ["SolvingSet","solution"])


text:
norm(A::Array; dims::Int) = sqrt.(sum(x->x^2, A, dims=dims))

function compute_r_nearest_image!(mdsim::MDSimulation)
	P = mdsim.params
	x = mdsim.state.x .% P.L

	mdsim.state.r_nearest_image = reshape(x', (1, P.d, P.N)) .- reshape(x, (P.N, P.d, 1))
	mdsim.state.r_nearest_image[mdsim.state.r_nearest_image .> P.L / 2]  .-= P.L
	mdsim.state.r_nearest_image[mdsim.state.r_nearest_image .< -P.L / 2] .+= P.L

	mdsim.state.r_nearest_image_magnitude = norm(mdsim.state.r_nearest_image, dims=2)
	mdsim.state.r_nearest_image_magnitude[P.diag_indices] .= Inf
end

function compute_row_vector_magnitudes(r)
	if size(r, 2) == 2
		return @. sqrt(r[:, 1]^2 + r[:, 2]^2)
	elseif size(r, 2) == 3
		return @. sqrt(r[:, 1]^2 + r[:, 2]^2 + r[:, 3]^2)
	end
end

function compute_energies!(mdsim::MDSimulation)
	P = mdsim.params
	α = (P.σ_LJ ./ mdsim.state.r_nearest_image_magnitude).^6
	mdsim.state.E_pot = sum((4 * P.ε_LJ) .* (α.^2 .- α)) / 2
	mdsim.state.E_kin = 0.5P.m * sum(mdsim.state.v.^2)
end

function compute_accelerations!(mdsim::MDSimulation)
	P = mdsim.params
	mdsim.state.a_prev = mdsim.state.a
	compute_r_nearest_image!(mdsim)

	α = (P.σ_LJ ./ mdsim.state.r_nearest_image_magnitude).^6
	mdsim.state.a = sum(
		((48 * P.ε_LJ / P.m) * (α.^2 - 0.5α) ./ mdsim.state.r_nearest_image_magnitude.^2)
		.* mdsim.state.r_nearest_image,
		dims=1
	)[1, :, :]'
end