text:
include("RouteguideClients/RouteguideClients.jl")
using .RouteguideClients

Base.show(io::IO, location::RouteguideClients.Point) = print(io, string("[", location.latitude, ", ", location.longitude, "]"))
Base.show(io::IO, feature::RouteguideClients.Feature) = print(io, string(feature.name, " - ", feature.location))
Base.show(io::IO, summary::RouteguideClients.RouteSummary) = print(io, string(summary.point_count, " points, ", summary.feature_count, " features, distance=", summary.distance, ", elapsed_time=", summary.elapsed_time))
Base.show(io::IO, note::RouteguideClients.RouteNote) = print(io, string(note.message, " ", note.location))

function randomPoint()
    latitude = (abs(rand(Int) % 180) - 90) * 1e7
    longitude = (abs(rand(Int) % 360) - 180) * 1e7
    RouteguideClients.Point(; latitude=latitude, longitude=longitude)
end

# single request, single response
function test_get_feature(client::RouteGuideBlockingClient)
    # existing feature
    point = RouteguideClients.Point(; latitude=409146138, longitude=-746188906)
    feature, status_future = RouteguideClients.GetFeature(client, point)
    gRPCCheck(status_future)
    @test isa(feature, RouteguideClients.Feature)
    @debug("existing feature", feature)

    # missing feature
    point = RouteguideClients.Point(; latitude=0, longitude=0)
    feature, status_future = RouteguideClients.GetFeature(client, point)
    gRPCCheck(status_future)
    @test isa(feature, RouteguideClients.Feature)
    @debug("missing feature", feature)
end

# single request, streaming response
function test_list_features(client::RouteGuideBlockingClient)
    @debug("listing features in an area")
    rect = RouteguideClients.Rectangle(; lo=RouteguideClients.Point(; latitude=400000000, longitude=-750000000), hi=RouteguideClients.Point(; latitude=420000000, longitude=-730000000))
    features, status_future = RouteguideClients.ListFeatures(client, rect)
    while isopen(features) || isready(features)
        try
            feature = take!(features)
            @debug(feature)
        catch ex
            (!isa(ex, InvalidStateException) || !fetch(status_future).success) && rethrow(ex)
        end
    end
    gRPCCheck(status_future)
    @test isa(features, Channel{RouteguideClients.Feature})
    @test !isopen(features)
end

# streaming request, single response
function test_record_route(client::RouteGuideBlockingClient)
    @sync begin
        point_count = abs(rand(Int) % 100) + 2
        @debug("recording a route", point_count)
        points_channel = Channel{RouteguideClients.Point}(1)
        @async begin
            for idx in 1:point_count
                put!(points_channel, randomPoint())
            end
            close(points_channel)
        end
        route_summary, status_future = RouteguideClients.RecordRoute(client, points_channel)
        gRPCCheck(status_future)
        @test isa(route_summary, RouteguideClients.RouteSummary)
        @test !isopen(points_channel)
        @debug("route summary: $route_summary")
    end
end

# streaming request, streaming response
function test_route_chat(client::RouteGuideBlockingClient)
    @sync begin
        notes = RouteguideClients.RouteNote[
            RouteguideClients.RouteNote(;location=RouteguideClients.Point(;latitude=0, longitude=1), message="First message"),
            RouteguideClients.RouteNote(;location=RouteguideClients.Point(;latitude=0, longitude=2), message="Second message"),
            RouteguideClients.RouteNote(;location=RouteguideClients.Point(;latitude=0, longitude=3), message="Third message"),
            RouteguideClients.RouteNote(;location=RouteguideClients.Point(;latitude=0, longitude=1), message="Fourth message"),
            RouteguideClients.RouteNote(;location=RouteguideClients.Point(;latitude=0, longitude=2), message="Fifth message"),
            RouteguideClients.RouteNote(;location=RouteguideClients.Point(;latitude=0, longitude=3), message="Sixth message"),
        ]
        @debug("route chat")
        in_channel = Channel{RouteguideClients.RouteNote}(1)
        @async begin
            for note in notes
                put!(in_channel, note)
            end
            close(in_channel)
        end
        out_channel, status_future = RouteguideClients.RouteChat(client, in_channel)
        nreceived = 0
        for note in out_channel
            nreceived += 1
            @debug("received note $note")
        end
        gRPCCheck(status_future)
        @test nreceived > 0
        @test isa(out_channel, Channel{RouteguideClients.RouteNote})
        @test !isopen(out_channel)
        @test !isopen(in_channel)
    end
end

function test_exception()
    client = RouteGuideBlockingClient("https://localhost:30000"; verbose=false)
    point = RouteguideClients.Point(; latitude=409146138, longitude=-746188906)
    @test_throws gRPCServiceCallException RouteguideClients.GetFeature(client, point)

    @test_throws ArgumentError RouteGuideBlockingClient("https://localhost:30000"; maxage=-1)
end

function test_message_length_limit(server_endpoint)
    point = RouteguideClients.Point(; latitude=409146138, longitude=-746188906)

    # test send message length limit
    client = RouteGuideBlockingClient(server_endpoint; max_message_length=1, verbose=false)
    @test_throws gRPCMessageTooLargeException RouteguideClients.GetFeature(client, point)

    # test recv message length limit
    client = RouteGuideBlockingClient(server_endpoint; max_recv_message_length=1, verbose=false)
    @test_throws gRPCMessageTooLargeException RouteguideClients.GetFeature(client, point)

    iob = IOBuffer()
    show(iob, gRPCMessageTooLargeException(1, 2))
    @test String(take!(iob)) == "gRPMessageTooLargeException(1, 2) - Encountered message size 2 > max configured 1"
    show(iob, gRPCServiceCallException(0, "test error"))
    @test String(take!(iob)) == "gRPCServiceCallException: 0, test error"
end

function test_async_get_feature(client::RouteGuideClient)
    # existing feature
    point = RouteguideClients.Point(; latitude=409146138, longitude=-746188906)
    results = Channel{Any}(1)
    RouteguideClients.GetFeature(client, point, result->put!(results, result))
    feature, status_future = take!(results)
    gRPCCheck(status_future)
    @test isa(feature, RouteguideClients.Feature)
    @debug("existing feature", feature)
end

function test_async_client(server_endpoint::String)
    client = RouteGuideClient(server_endpoint; verbose=false)
    @testset "GetFeature" begin
        test_async_get_feature(client)
    end
end

function test_blocking_client(server_endpoint::String)
    client = RouteGuideBlockingClient(server_endpoint; verbose=false)
    @testset "request response" begin
        test_get_feature(client)
    end
    @testset "streaming recv" begin
        test_list_features(client)
    end
    @testset "streaming send" begin
        test_record_route(client)
    end
    @testset "streaming send recv" begin
        test_route_chat(client)
    end
    @testset "error handling" begin
        test_exception()
    end
    @testset "message length limits" begin
        test_message_length_limit(server_endpoint)
    end
end

function test_clients(server_endpoint::String)
    @info("testing blocking client")
    test_blocking_client(server_endpoint)
    @info("testing async client")
    test_async_client(server_endpoint)
end

function test_task_safety(server_endpoint::String)
    @testset "async safety" begin
        client = RouteGuideBlockingClient(server_endpoint; verbose=false)
        @sync begin
            for taskid in 1:5
                @async begin
                    test_get_feature(client)
                    test_list_features(client)
                    test_record_route(client)
                    test_route_chat(client)
                    test_exception()
                    test_message_length_limit(server_endpoint)
                end
            end
        end
    end
end


text:
function Shader(f::File{format"GLSLShader"})
    st = stream(open(f))
    s = Shader(Symbol(f.filename), read(st), shadertype(f))
    close(st)
    s
end
Shader(s::Shader; name=s.name, source=s.source, typ=s.typ) = Shader(name, source, typ)
# Different shader string literals- usage: e.g. frag" my shader code"
macro frag_str(source::AbstractString)
    quote
        Shader(Symbol(@__FILE__), $(Vector{UInt8}(ascii(source))), GL_FRAGMENT_SHADER)
    end
end
macro vert_str(source::AbstractString)
    quote
        Shader(Symbol(@__FILE__), $(Vector{UInt8}(ascii(source))), GL_VERTEX_SHADER)
    end
end
macro geom_str(source::AbstractString)
    quote
        Shader(Symbol(@__FILE__), $(Vector{UInt8}(ascii(source))), GL_GEOMETRY_SHADER)
    end
end
macro comp_str(source::AbstractString)
    quote
        Shader(Symbol(@__FILE__), $(Vector{UInt8}(ascii(source))), GL_COMPUTE_SHADER)
    end
end

function getinfolog(obj::GLuint)
    # Return the info log for obj, whether it be a shader or a program.
    isShader    = glIsShader(obj)
    getiv       = isShader == GL_TRUE ? glGetShaderiv : glGetProgramiv
    get_log     = isShader == GL_TRUE ? glGetShaderInfoLog : glGetProgramInfoLog

    # Get the maximum possible length for the descriptive error message
    maxlength = GLint[0]
    getiv(obj, GL_INFO_LOG_LENGTH, maxlength)
    maxlength = first(maxlength)
    # Return the text of the message if there is any
    if maxlength > 0
        buffer = zeros(GLchar, maxlength)
        sizei = GLsizei[0]
        get_log(obj, maxlength, sizei, buffer)
        length = first(sizei)
        return unsafe_string(pointer(buffer), length)
    else
        return "success"
    end
end

function iscompiled(shader::GLuint)
    success = GLint[0]
    glGetShaderiv(shader, GL_COMPILE_STATUS, success)
    return first(success) == GL_TRUE
end
islinked(program::GLuint) = glGetProgramiv(program, GL_LINK_STATUS) == GL_TRUE

function createshader(shadertype::GLenum)
    shaderid = glCreateShader(shadertype)
    @assert shaderid > 0 "opengl context is not active or shader type not accepted. Shadertype: $(GLENUM(shadertype).name)"
    shaderid::GLuint
end
function createprogram()
    program = glCreateProgram()
    @assert program > 0 "couldn't create program. Most likely, opengl context is not active"
    program::GLuint
end

shadertype(s::Shader)                      = s.typ
function shadertype(f::File{format"GLSLShader"})
    ext = file_extension(f)
    ext == ".comp" && return GL_COMPUTE_SHADER
    ext == ".vert" && return GL_VERTEX_SHADER
    ext == ".frag" && return GL_FRAGMENT_SHADER
    ext == ".geom" && return GL_GEOMETRY_SHADER
    error("$ext not a valid extension for $f")
end

"""
Gives back a signal, which signals true everytime the file gets edited
"""
function isupdated(file::File, updatewhile=Signal(true), update_interval=1.0)
    fn = filename(file)
    file_edited = foldp((false, mtime(fn)), fpswhen(updatewhile, 1.0/update_interval)) do v0, v1
        time_edited = mtime(fn)
        (!isapprox(0.0, v0[2] - time_edited), time_edited)
    end
    preserve(file_edited)
    return preserve(filter(identity, false, const_lift(first, file_edited))) # extract bool
end

#reads from the file and updates the source whenever the file gets edited
function const_lift_shader(shader_file::File, updatewhile=Signal(true), update_interval=1.0)
    s = const_lift(isupdated(shader_file, updatewhile, update_interval)) do _unused
        Shader(shader_file)
    end
    preserve(s)
end

#Implement File IO interface
load(f::File{format"GLSLShader"}) = preserve(const_lift_shader(f))
function save(f::File{format"GLSLShader"}, data::Shader)
    s = open(f, "w")
    write(s, data.source)
    close(s)
end

compileshader(file::File{format"GLSLShader"}, program::GLuint) = compileshader(load(file), program)
                    #(shadertype, shadercode) -> shader id
let shader_cache = Dict{Tuple{GLenum, Vector{UInt8}}, GLuint}() # shader cache prevents that a shader is compiled more than one time
    #finalizer(shader_cache, dict->foreach(glDeleteShader, values(dict))) # delete all shaders when done
    empty_shader_cache!() = empty!(shader_cache)
    global empty_shader_cache!

    function compileshader(shader::Shader)
        get!(shader_cache, (shader.typ, shader.source)) do
            shaderid = createshader(shader.typ)
            glShaderSource(shaderid, shader.source)
            glCompileShader(shaderid)
            if !iscompiled(shaderid)
                print_with_lines(Compat.String(shader.source))
                warn("shader $(shader.name) didn't compile. \n$(getinfolog(shaderid))")
            end
            shaderid
        end
    end
end

export empty_shadercache

function uniformlocations(nametypedict::Dict{Symbol, GLenum}, program)
    isempty(nametypedict) && return Dict{Symbol,Tuple}()
    texturetarget = -1 # start -1, as texture samplers start at 0
    return Dict{Symbol, Tuple}(map(nametypedict) do name_type
        name, typ = name_type
        loc = get_uniform_location(program, name)
        str_name = string(name)
        if istexturesampler(typ)
            texturetarget += 1
            return name => (loc, texturetarget)
        else
            return name => (loc,)
        end
    end)
end

# Actually compiles and links shader sources
function GLProgram(
        shaders::Vector{Shader}, program=createprogram();
        fragdatalocation=Tuple{Int, Compat.UTF8String}[]
    )

    # Remove old shaders
    glUseProgram(0)
    shader_ids = glGetAttachedShaders(program)
    foreach(glDetachShader, program, shader_ids)

    #attach new ones
    shader_ids = map(shaders) do shader
        shaderid = compileshader(shader)
        glAttachShader(program, shaderid)
        shaderid
    end

    #Bind frag data
    for (location, name) in fragdatalocation
        glBindFragDataLocation(program, location, ascii(name))
    end

    #link program
    glLinkProgram(program)
    !islinked(program) && warn("program $program not linked. Error in: \n", join(map(x->x.name, shaders), " or\n"), "\n", getinfolog(program))
    #foreach(glDeleteShader, shader_ids) # Can be deleted, as they will still be linked to Program and released after program gets released

    # generate the link locations
    nametypedict        = uniform_name_type(program)
    uniformlocationdict = uniformlocations(nametypedict, program)

    GLProgram(program, shaders, nametypedict, uniformlocationdict)
end

abstract AbstractLazyShader
immutable LazyShader <: AbstractLazyShader
    paths  ::Tuple
    kw_args::Vector
    function LazyShader(paths...; kw_args...)
        new(paths, kw_args)
    end
end

gl_convert(lazyshader::AbstractLazyShader, data) = TemplateProgram(
    lazyshader.paths...;
    attributes = data,
    lazyshader.kw_args...
)

function insert_from_view(io, replace_view::Function, keyword::AbstractString)
    print(io, replace_view(keyword))
    nothing
end

function insert_from_view(io, replace_view::Dict, keyword::AbstractString)
    if haskey(replace_view, keyword)
        print(io, replace_view[keyword])
    end
    nothing
end

"""
Replaces
{{keyword}} with the key in `replace_view`, or replace_view(key)
in a string
"""
function mustache_replace(replace_view::Union{Dict, Function}, string)
    io = IOBuffer()
    replace_started = false
    open_mustaches = 0
    closed_mustaches = 0
    i = 0
    replace_begin = i
    last_char = SubString(string, 1, 1)
    len = endof(string)
    while i <= len
        i = nextind(string, i)
        char = SubString(string, i, i)
        if replace_started
            # ignore, or wait for }
            if char == "}"
                closed_mustaches += 1
                if closed_mustaches == 2 # we found a complete mustache!
                    insert_from_view(io, replace_view, SubString(string, replace_begin+1, i-2))
                    open_mustaches = 0
                    closed_mustaches = 0
                    replace_started = false
                end
            else
                closed_mustaches = 0
                continue
            end
        elseif char == "{"
            open_mustaches += 1
            if open_mustaches == 2
                replace_begin = i
                replace_started = true
            end
        else
            if open_mustaches == 1
                print(io, last_char)
            end
            print(io, char) # just copy all the rest
            open_mustaches = 0
            closed_mustaches = 0
        end
        last_char = char
    end
    takebuf_string(io)
end

function mustache2replacement(mustache_key, view, attributes)
    haskey(view, mustache_key) && return view[mustache_key]
    for postfix in ("_type", "_calculation")
        keystring = replace(mustache_key, postfix, "")
        keysym = Symbol(keystring)
        if haskey(attributes, keysym)
            val = attributes[keysym]
            if !isa(val, AbstractString)
                return if postfix == "_type"
                    toglsltype_string(val)
                else  postfix == "_calculation"
                    glsl_variable_access(keystring, val)
                end
            end
        end
    end
    "" # no match found, leave empty!
end
# Takes a shader template and renders the template and returns shader source
template2source(source::Array{UInt8, 1}, attributes::Dict{Symbol, Any}, view) = template2source(Compat.String(source), attributes, view)
function template2source(source::AbstractString, attributes::Dict{Symbol, Any}, view)
    source = mustache_replace(source) do mustache_key
        mustache2replacement(mustache_key, view, attributes)
    end
    source
end

#TemplateProgram() = error("Can't create TemplateProgram without parameters")

function TemplateProgram(x::Union{Shader, File, Signal{Shader}}...; kw_args...)
    TemplateProgram(merge(Dict(
        :view               => Dict{Compat.UTF8String, Compat.UTF8String}(),
        :attributes         => Dict{Symbol, Any}(),
        :fragdatalocation   => Tuple{Int, Compat.UTF8String}[],
        :program            => createprogram()
    ), Dict{Symbol, Any}(kw_args)), x...)
end

function TemplateProgram(kw_args::Dict{Symbol, Any}, s::File, shaders::File...)
    updatewhile     = get(kw_args, :updatewhile, Signal(true))
    update_interval = get(kw_args, :update_interval, 1.0)
    shader_signals  = map(s->const_lift_shader(s, updatewhile, update_interval), [s,shaders...])
    TemplateProgram(kw_args, shader_signals...)
end
function TemplateProgram(kw_args::Dict{Symbol, Any}, s::Signal{Shader}, shaders::Signal{Shader}...)
    program_signal = const_lift(s, shaders...) do _unused... #just needed to update the signal
        # extract values from signals
        shader_values = map(value, [s, shaders...])
        TemplateProgram(kw_args, shader_values...)
    end
    Reactive.preserve(program_signal)
    program_signal
end
function TemplateProgram(kw_args::Dict{Symbol, Any}, shaders::Union{Shader,Signal{Shader}}...)
    newshaders = map(s->isa(s, Shader) ? Signal(s) : s, shaders)
    TemplateProgram(kw_args, newshaders...)
end

function TemplateProgram(kw_args::Dict{Symbol, Any}, s::Shader, shaders::Shader...)
    @materialize program, view, attributes, fragdatalocation = kw_args
    if haskey(view, "in") || haskey(view, "out") || haskey(view, "GLSL_VERSION")
        println("warning: using internal keyword \"$(in/out/GLSL_VERSION)\" for shader template. The value will be overwritten")
    end
    extension = is_apple() ? "" : "#extension GL_ARB_draw_instanced : enable"
    if haskey(view, "GLSL_EXTENSIONS")
        #to do: check custom extension...
        #for now we just append the extensions
        extension *= "\n" * view["GLSL_EXTENSIONS"]
    end
    internaldata = Dict{Compat.UTF8String, Compat.UTF8String}(
        "GLSL_VERSION"    => glsl_version_string(),
        "GLSL_EXTENSIONS" => extension
    )
    view = merge(internaldata, view)

    # transform dict of templates into actual shader source
    code = Shader[Shader(shader, source=template2source(shader.source, attributes, view)) for shader in [s, shaders...]]
    return GLProgram(code, program, fragdatalocation=fragdatalocation)
end

function glsl_version_string()
    glsl = split(unsafe_string(glGetString(GL_SHADING_LANGUAGE_VERSION)), ['.', ' '])
    if length(glsl) >= 2
        glsl = VersionNumber(parse(Int, glsl[1]), parse(Int, glsl[2]))
        glsl.major == 1 && glsl.minor <= 2 && error("OpenGL shading Language version too low. Try updating graphic driver!")
        glsl_version = string(glsl.major) * rpad(string(glsl.minor),2,"0")
        return "#version $(glsl_version)\n"
    else
        error("could not parse GLSL version: $glsl")
    end
end

ext: jl


text:
module CircularArrayBuffers

export CircularArrayBuffer, CircularVectorBuffer, capacity, isfull

"""
    CircularArrayBuffer{T}(sz::Integer...) -> CircularArrayBuffer{T, N}

`CircularArrayBuffer` uses a `N`-dimension `Array` of size `sz` to serve as a buffer for
`N-1`-dimension `Array`s of the same size.
"""
mutable struct CircularArrayBuffer{T,N} <: AbstractArray{T,N}
    buffer::Array{T,N}
    first::Int
    nframes::Int
    step_size::Int
end

const CircularVectorBuffer{T} = CircularArrayBuffer{T,1}

CircularVectorBuffer{T}(n::Integer) where {T} = CircularArrayBuffer{T}(n)

function CircularArrayBuffer{T}(d::Integer...) where {T}
    N = length(d)
    CircularArrayBuffer(Array{T}(undef, d...), 1, 0, N == 1 ? 1 : *(d[1:end-1]...))
end

function CircularArrayBuffer(A::AbstractArray{T,N}) where {T,N}
    CircularArrayBuffer(A, 1, size(A, N), N == 1 ? 1 : *(size(A)[1:end-1]...))
end

Base.IndexStyle(::CircularArrayBuffer) = IndexLinear()

Base.size(cb::CircularArrayBuffer{T,N}, i::Integer) where {T,N} = i == N ? cb.nframes : size(cb.buffer, i)
Base.size(cb::CircularArrayBuffer{T,N}) where {T,N} = ntuple(i -> size(cb, i), N)
Base.getindex(cb::CircularArrayBuffer{T,N}, i::Int) where {T,N} = getindex(cb.buffer, _buffer_index(cb, i))
Base.setindex!(cb::CircularArrayBuffer{T,N}, v, i::Int) where {T,N} = setindex!(cb.buffer, v, _buffer_index(cb, i))

capacity(cb::CircularArrayBuffer{T,N}) where {T,N} = size(cb.buffer, N)
isfull(cb::CircularArrayBuffer) = cb.nframes == capacity(cb)
Base.isempty(cb::CircularArrayBuffer) = cb.nframes == 0

@inline function _buffer_index(cb::CircularArrayBuffer, i::Int)
    ind = (cb.first - 1) * cb.step_size + i
    if ind > length(cb.buffer)
        ind - length(cb.buffer)
    else
        ind
    end
end

@inline function _buffer_frame(cb::CircularArrayBuffer, i::Int)
    n = capacity(cb)
    idx = cb.first + i - 1
    if idx > n
        idx - n
    else
        idx
    end
end

_buffer_frame(cb::CircularArrayBuffer, I::Vector{Int}) = map(i -> _buffer_frame(cb, i), I)

function Base.empty!(cb::CircularArrayBuffer)
    cb.nframes = 0
    cb
end

function Base.push!(cb::CircularArrayBuffer{T,N}, data) where {T,N}
    if cb.nframes == capacity(cb)
        cb.first = (cb.first == capacity(cb) ? 1 : cb.first + 1)
    else
        cb.nframes += 1
    end
    if N == 1
        if ndims(data) == 0
            cb[cb.nframes] = data[]
        else
            cb[cb.nframes] = data
        end
    else
        cb[ntuple(_ -> (:), N - 1)..., cb.nframes] .= data
    end
    cb
end

function Base.append!(cb::CircularArrayBuffer{T,N}, data) where {T,N}
    d, r = divrem(length(data), cb.step_size)
    @assert r == 0
    if length(data) >= length(cb.buffer)
        cb.nframes = capacity(cb)
        cb.first = 1
        cb.buffer[:] .= @view data[end-length(cb.buffer)+1:end]
    else
        start_idx = (cb.first - 1) * cb.step_size + length(cb) + 1
        end_idx = start_idx + length(data) - 1
        if start_idx > length(cb.buffer)
            start_idx -= length(cb.buffer)
            end_idx -= length(cb.buffer)
        end
        if end_idx > length(cb.buffer)
            n_first_part = length(cb.buffer) - start_idx + 1
            n_second_part = length(data) - n_first_part
            cb.buffer[end-n_first_part+1:end] .= @view data[1:n_first_part]
            cb.buffer[1:n_second_part] .= @view data[end-n_second_part+1:end]
        else
            cb.buffer[start_idx:end_idx] .= data
        end

        if cb.nframes + d > capacity(cb)
            cb.first += cb.nframes + d - capacity(cb)
            if cb.first > capacity(cb)
                cb.first -= capacity(cb)
            end
            cb.nframes = capacity(cb)
        else
            cb.nframes += d
        end
    end
    cb
end

function Base.pop!(cb::CircularArrayBuffer{T,N}) where {T,N}
    if cb.nframes <= 0
        throw(ArgumentError("buffer must be non-empty"))
    else
        res = @views cb.buffer[ntuple(_ -> (:), N - 1)..., _buffer_frame(cb, cb.nframes)]
        cb.nframes -= 1
        res
    end
end

function Base.popfirst!(cb::CircularArrayBuffer{T,N}) where {T,N}
    if cb.nframes <= 0
        throw(ArgumentError("buffer must be non-empty"))
    else
        res = @views cb.buffer[ntuple(_ -> (:), N - 1)..., _buffer_frame(cb, 1)]
        cb.nframes -= 1
        cb.first += 1
        if cb.first > capacity(cb)
            cb.first = 1
        end
        res
    end
end

end


text:
@testset "Jordan and Rational canonical form" begin

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^6 + rand(Qx, 1:5, -3:3)
  M = companion_matrix(f)
  @test minpoly(Qx, M) == f
  @test charpoly(Qx, M) == f

  ff = factor(f)
  g = first(keys(ff.fac))
  g = divexact(g, leading_coefficient(g))
  J = Hecke.jordan_block(g, 2)
  @test minpoly(Qx, J) == g^2
  @test charpoly(Qx, J) == g^2

  M = identity_matrix(FlintQQ, 6)
  for i = 1:6
    for j = i+1:6
      M[i, j] = rand(FlintQQ, -2:2)
    end
  end
  M1 = transpose(M)
  @test Hecke.is_similar(M, M1)
  MT = Hecke.conjugating_matrix(M, M1)
  @test MT*M1*inv(MT) == M
  J, S = jordan_normal_form(M)
  @test S*M*inv(S) == J
  for i = 1:6
    for j = 1:i-1
      @test iszero(J[i, j])
    end
  end
  for i = 1:5
    @test isone(J[i, i+1]) || iszero(J[i, i+1])
  end
  for i = 1:6
    for j = i+2:6
      @test iszero(J[i, i+2])
    end
  end

  f = x^3+ rand(Qx, 1:2, -2:2)
  g = x^3 + rand(Qx, 1:2, -2:2)
  f1 = f
  f2 = f*g
  C = zero_matrix(FlintQQ, 9, 9)
  Hecke._copy_matrix_into_matrix(C, 1, 1, companion_matrix(f1))
  Hecke._copy_matrix_into_matrix(C, 4, 4, companion_matrix(f2))
  CF, TM = rational_canonical_form(C)
  @test CF == C
  @test TM * C * inv(TM) == CF

end

@testset "Spectrum and eigenspaces" begin
  M = matrix(FlintQQ, 4, 4, [ 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 ])
  l = eigvals(M)
  @test length(keys(l)) == 1
  @test haskey(l, one(FlintQQ))
  @test l[one(FlintQQ)] == 2

  K, a = CyclotomicField(3, "a")
  lK = eigvals(M, K)
  @test length(keys(lK)) == 3
  @test haskey(lK, one(K)) && haskey(lK, a) && haskey(lK, -a - 1)
  @test lK[one(K)] == 2
  @test lK[a] == 1
  @test lK[-a - 1] == 1

  M = matrix(K, 4, 4, [ 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 ])
  lK2 = eigvals(M)
  @test lK == lK2

  Eig = eigenspaces(M, side = :right)
  V = zero_matrix(K, 4, 0)
  for (e, v) in Eig
    @test haskey(lK2, e)
    @test lK2[e] == ncols(v)
    @test M*v == e*v
    V = hcat(V, v)
  end
  @test rref!(V) == 4

  Eig = eigenspaces(M, side = :left)
  V = zero_matrix(K, 0, 4)
  for (e, v) in Eig
    @test haskey(lK2, e)
    @test lK2[e] == nrows(v)
    @test v*M == e*v
    V = vcat(V, v)
  end
  @test rref!(V) == 4

  M = [ matrix(QQ, [ 0 1 0 0 0 0;
                    1 0 0 0 0 0;
                    0 0 1 0 0 0;
                    0 0 0 1 0 0;
                    0 0 0 0 1 0;
                    0 0 0 0 0 1; ]),
       matrix(QQ, [ 1 0 0 0 0 0;
                    0 1 0 0 0 0;
                    0 0 0 1 0 0;
                    0 0 1 0 0 0;
                    0 0 0 0 1 0;
                    0 0 0 0 0 1; ]),
       matrix(QQ, [ 1 0 0 0 0 0;
                    0 1 0 0 0 0;
                    0 0 1 0 0 0;
                    0 0 0 1 0 0;
                    0 0 0 0 0 1;
                    0 0 0 0 1 0; ]) ]

  Eig = common_eigenspaces(M, side = :right)
  V = zero_matrix(FlintQQ, 6, 0)
  for (e, v) in Eig
    @test length(e) == 3
    for i = 1:3
      @test M[i]*v == e[i]*v
    end
    V = hcat(V, v)
  end
  @test rref!(V) == 6

  Eig = common_eigenspaces(M, side = :left)
  V = zero_matrix(K, 0, 6)
  for (e, v) in Eig
    @test length(e) == 3
    for i = 1:3
      @test v*M[i] == e[i]*v
    end
    V = vcat(V, v)
  end
  @test rref!(V) == 6

  N = [ matrix(QQ, [ 0 0 1 0 0 0;
                     1 0 0 0 0 0;
                     0 1 0 0 0 0;
                     0 0 0 1 0 0;
                     0 0 0 0 1 0;
                     0 0 0 0 0 1; ]),
        matrix(QQ, [ 1 0 0 0 0 0;
                     0 1 0 0 0 0;
                     0 0 1 0 0 0;
                     0 0 0 0 0 1;
                     0 0 0 1 0 0;
                     0 0 0 0 1 0; ]) ]

  @test_throws ErrorException Hecke.simultaneous_diagonalization(N)

  T, D = Hecke.simultaneous_diagonalization(N, K)
  for i = 1:length(N)
    @test T*N[i]*inv(T) == D[i]
    @test is_invertible(D[i])
    @test is_diagonal(D[i])
  end
end


text:
# A bunch of functions to perform a single iteration to minimize h(t)
# and/or compute a root of h'(t)
include("DividedDifference2.jl")
include("coefHIN.jl")

# basic iterations
base =  Function[]

"""One iteration of the secant method"""
function s_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    s = xm - x
    y = dhm - dh

    d = - (dh * s / y)

    xp = x + d
    return xp
end
push!(base, s_iteration)

"""One iteration of the Newton's method"""
function N_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    d = - (dh / d2h)

    xp = x + d
    return xp
end
push!(base, N_iteration)

"""One iteration of a cubic Taylor approximation
   may also be described as a quadratic Taylor approximation of h'
"""
function C_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,d3h,d3hm)
    # solution cubique...
    disc = d2h^2 - 2.0*d3h*dh
    d = disc < 0 ? -sign(d3h)*Inf : (-d2h + sqrt(disc)) / d3h

    xp = x + d
    return xp
end
push!(base, C_iteration)

# iterations on the quintic interpolant of h, h' and h''
quintic = Function[]

"""Two Newton iteration on the quintic interpolant of h
   at the current and previous point
"""
function N2Opt_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    d = - (dh / d2h)
    x2 = x + d
    # estimate dh2 from the hermite interpolant

    a = coefHIN5([x;xm],[h dh d2h; hm dhm d2hm])

    p = Polynomial(a)
    dp = Polynomials.derivative(p)
    d2p = Polynomials.derivative(dp)

    d2 = -(dp(d) / d2p(d))

    xp = x2 + d2

    return xp
end
push!(quintic, N2Opt_iteration)

"""Three Newton iteration on the quintic interpolant of h
   at the current and previous point
"""
function N3Opt_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    d = - (dh / d2h)
    x2 = x + d
    # estimate dh2 from the hermite interpolant

    a = coefHIN5([x;xm],[h dh d2h; hm dhm d2hm])

    p = Polynomial(a)
    dp = Polynomials.derivative(p)
    d2p = Polynomials.derivative(dp)

    d2 = -(dp(d) / d2p(d))

    x3 = x2 + d2
    d3 = -(dp(d+d2) / d2p(d+d2))

    xp = x3 + d3
    return xp
end
push!(quintic, N3Opt_iteration)

# iterations on the cubic interpolant of h' and h''
cubic1 = Function[]

"""Two Newton iteration on the cubic interpolant of h'
   at the current and previous point
"""
function N2_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)

    d = - (dh / d2h)
    x2 = x + d
    # estimate dh2 from the hermite interpolant

    a = coefHIN3([x;xm],[dh d2h; dhm d2hm])
    p = Polynomial(a)
    dp = Polynomials.derivative(p)

    d2 = -(p(d) / dp(d))

    xp = x2 + d2
    return xp
end
push!(cubic1, N2_iteration)

"""Three Newton iterations on the cubic interpolant of h'
   at the current and previous point
"""
function N3_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    d = - (dh / d2h)
    x2 = x + d
    # estimate dh2 from the hermite interpolant

    a = coefHIN3([x;xm],[dh d2h; dhm d2hm])

    p = Polynomial(a)
    dp = Polynomials.derivative(p)

    d2 = -(p(d) / dp(d))

    x3 = x2 + d2
    d3 = -(p(d+d2) / dp(d+d2))
    xp = x3 + d3

    return xp
end
push!(cubic1, N3_iteration)

"""Two secant iterations on the cubic interpolant of h'
   at the current and previous point
"""
function s2Opt_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    s = x - xm
    y = dh - dhm

    d = - (dh * s / y)

    x2 = x + d
    # estimate dh2 from the hermite interpolant
    a = coefHIN3([x;xm],[h dh; hm dhm])
    s2 = x2 - x

    p = Polynomial(a)
    dp = Polynomials.derivative(p)
    dh2 = dp( s2)

    # perform a secant iteration with this estimated dh2
    y2 = dh2 - dh
    d2 = - (dh2 * s2 / y2)

    xp = x2 + d2
    if ( xp == x ) xp = x2 end #safeguard
    return xp
end
push!(cubic1, s2Opt_iteration)

# iterations on the cubic interpolant of h and h'
cubic0 = Function[]

"""Minimizer of the cubic interpolant of h and dh"""
function sCOpt_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    # Hermite cubic interpolant
    a = coefHIN3([x;xm],[h dh; hm dhm])
    p = Polynomial(a)
    dp = Polynomials.derivative(p)
    dh  = a[2]
    d2h = 2.0*a[3]
    d3h = 3.0*a[4]
    disc = d2h^2 - 2.0*d3h*dh
    d = disc < 0 ? -sign(d3h)*Inf : (-d2h + sqrt(disc)) / d3h

    xp = x + d

    return xp
end
push!(cubic0, sCOpt_iteration)

"""Two Newton's iterations on the cubic interpolant of h
   at the current and previous point
"""
function sNOpt_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    # Hermite cubic interpolant
    a = coefHIN3([x;xm],[h dh; hm dhm])
    p = Polynomial(a)
    dp = Polynomials.derivative(p)
    d2p = Polynomials.derivative(dp)

    # Newton step on the interpolant
    # Do not forget that the interpolant is rooted at 0
    d = - dp(0.0) / d2p(0.0)
    d2 = - dp(d) / d2p(d)
    xp = x + d + d2

    return xp
end
push!(cubic0, sNOpt_iteration)

"""Three secant iterations on the cubic interpolant of h
   at the current and previous point; not efficient????
"""
function s3Opt_iteration(x,xm,h,hm,dh,dhm,d2h,d2hm,::Any,::Any)
    s = x - xm
    y = dh - dhm

    d = - (dh * s / y)

    x2 = x + d

    # estimate dh2 from the hermite interpolant
    a = coefHIN3([x;xm],[h dh; hm dhm])
    p = Polynomial(a)
    dp = Polynomials.derivative(p)

    s2 = x2 - x
    dh2 = dp(s2)
    # perform a secant iteration with this estimated dh2
    y2 = dh2 - dh
    d2 = - (dh2 * s2 / y2)
    x3 = x2 + d2

    s3 = x3 - x2
    dh3 = dp(s3)
    # perform a secant iteration with this estimated dh3
    y3 = dh3 - dh2
    d3 = - (dh3 * s3 / y3)

    xp = x3 + d3
    # not good, seems to fall back close to regular secant iteration
    if ( xp == x ) xp = x3 end#(x+xm)/2.0 end#@show a end
    return xp
end
push!(cubic0, s3Opt_iteration)

export base, quintic, cubic1, cubic0

for n in nameof.(base ∪ quintic ∪ cubic1 ∪ cubic0 )
    @eval export $n
end


text:
using OnlineNMF
using OnlineNMF: readcsv, writecsv
using Test
using Pkg
using DelimitedFiles
using Statistics
using Distributions

tmp = mktempdir()
julia = joinpath(Sys.BINDIR, "julia")
bindir = joinpath(dirname(pathof(OnlineNMF)), "..", "bin")

function testfilesize(remove::Bool, x...)
  for n = 1:length(x)
    @test filesize(x[n]) != 0
    if remove
      rm(x[n])
    end
  end
end

#####################################
X = Int64.(ceil.(rand(NegativeBinomial(1, 0.5), 300, 30)))
X[1:50, 1:10] .= 100*X[1:50, 1:10]
X[51:100, 11:20] .= 100*X[51:100, 11:20]
X[101:150, 21:30] .= 100*X[101:150, 21:30]
writecsv(joinpath(tmp, "X.csv"), X)

Y = Int64.(ceil.(rand(NegativeBinomial(1, 0.5), 300, 40)))
Y[1:50, 1:15] .= 100*Y[1:50, 1:15]
Y[51:100, 16:24] .= 100*Y[51:100, 16:24]
Y[101:150, 25:40] .= 100*Y[101:150, 25:40]
writecsv(joinpath(tmp, "Y.csv"), Y)

Z = Int64.(ceil.(rand(NegativeBinomial(1, 0.5), 300, 50)))
Z[1:50, 1:20] .= 100*Z[1:50, 1:20]
Z[51:100, 21:40] .= 100*Z[51:100, 21:40]
Z[101:150, 41:50] .= 100*Z[101:150, 41:50]
writecsv(joinpath(tmp, "Z.csv"), Z)
#####################################

#####################################
println("####### Binarization (Julia API) #######")
csv2bin(csvfile=joinpath(tmp, "X.csv"),
  binfile=joinpath(tmp, "X.zst"))
csv2bin(csvfile=joinpath(tmp, "Y.csv"),
  binfile=joinpath(tmp, "Y.zst"))
csv2bin(csvfile=joinpath(tmp, "Z.csv"),
  binfile=joinpath(tmp, "Z.zst"))

testfilesize(false, joinpath(tmp, "X.zst"))
testfilesize(false, joinpath(tmp, "Y.zst"))
testfilesize(false, joinpath(tmp, "Z.zst"))
####################################

# ####################################
# println("####### Binarization (Command line) #######")
# run(`$(julia) $(joinpath(bindir, "csv2bin")) --csvfile $(joinpath(tmp, "X.csv")) --binfile $(joinpath(tmp, "X.zst"))`)
# run(`$(julia) $(joinpath(bindir, "csv2bin")) --csvfile $(joinpath(tmp, "Y.csv")) --binfile $(joinpath(tmp, "Y.zst"))`)
# run(`$(julia) $(joinpath(bindir, "csv2bin")) --csvfile $(joinpath(tmp, "Z.csv")) --binfile $(joinpath(tmp, "Z.zst"))`)

# testfilesize(false, joinpath(tmp, "X.zst"))
# testfilesize(false, joinpath(tmp, "Y.zst"))
# testfilesize(false, joinpath(tmp, "Z.zst"))
# ####################################

# ####################################
# println("####### Multiplicative Update (Julia API) #######")
# # Alpha = 2 : Pearson divergence-based
# out_mu_alpha1 = mu(input=joinpath(tmp, "X.zst"), alpha=2,
#   dim=3, numepoch=10, numbatch=10, logdir=tmp)

# @test size(out_mu_alpha1[1][1]) == (30, 3)
# @test size(out_mu_alpha1[1][2]) == (40, 3)
# @test size(out_mu_alpha1[1][3]) == (50, 3)
# @test size(out_mu_alpha1[2][1]) == (3, )
# @test size(out_mu_alpha1[2][2]) == (3, )
# @test size(out_mu_alpha1[2][3]) == (3, )
# @test size(out_mu_alpha1[3][1]) == (300, 3)
# @test size(out_mu_alpha1[3][2]) == (300, 3)
# @test size(out_mu_alpha1[3][3]) == (300, 3)
# @test size(out_mu_alpha1[4][1]) == ()
# @test size(out_mu_alpha1[4][2]) == ()
# @test size(out_mu_alpha1[4][3]) == ()
# @test size(out_mu_alpha1[5][1]) == ()
# @test size(out_mu_alpha1[5][2]) == ()
# @test size(out_mu_alpha1[5][3]) == ()
# @test size(out_mu_alpha1[6][1]) == ()
# @test size(out_mu_alpha1[6][2]) == ()
# @test size(out_mu_alpha1[6][3]) == ()
# @test size(out_mu_alpha1[7][1]) == ()
# @test size(out_mu_alpha1[7][2]) == ()
# @test size(out_mu_alpha1[7][3]) == ()
# @test size(out_mu_alpha1[8]) == ()

# # Alpha = 1 : KL divergence-based

# # Alpha = 0.5 : Hellinger divergence-based

# # Alpha = -1 : Neyman divergence-based

# # Beta = 2 : Euclid distance-based

# # Beta = 0 : Itakura-Saito divergence-based

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

# ####################################
# println("####### Multiplicative Update (Command line) #######")

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

# ####################################
# println("####### Simultaneous NMF (Julia API) #######")
# # Alpha = 2 : Pearson divergence-based
# out_sinmf1 = sinmf(
#   input=[joinpath(tmp, "X.zst"),
#     joinpath(tmp, "Y.zst"),
#     joinpath(tmp, "Z.zst")],
#   alpha = 2,
#   dim=3,
#   numepoch=10,
#   numbatch=10,
#   logdir=tmp)

# @test size(out_sinmf1[1][1]) == (30, 3)
# @test size(out_sinmf1[1][2]) == (40, 3)
# @test size(out_sinmf1[1][3]) == (50, 3)
# @test size(out_sinmf1[2][1]) == (3, )
# @test size(out_sinmf1[2][2]) == (3, )
# @test size(out_sinmf1[2][3]) == (3, )
# @test size(out_sinmf1[3][1]) == (300, 3)
# @test size(out_sinmf1[3][2]) == (300, 3)
# @test size(out_sinmf1[3][3]) == (300, 3)
# @test size(out_sinmf1[4][1]) == ()
# @test size(out_sinmf1[4][2]) == ()
# @test size(out_sinmf1[4][3]) == ()
# @test size(out_sinmf1[5][1]) == ()
# @test size(out_sinmf1[5][2]) == ()
# @test size(out_sinmf1[5][3]) == ()
# @test size(out_sinmf1[6][1]) == ()
# @test size(out_sinmf1[6][2]) == ()
# @test size(out_sinmf1[6][3]) == ()
# @test size(out_sinmf1[7][1]) == ()
# @test size(out_sinmf1[7][2]) == ()
# @test size(out_sinmf1[7][3]) == ()
# @test size(out_sinmf1[8]) == ()

# # Alpha = 1 : KL divergence-based

# # Alpha = 0.5 : Hellinger divergence-based

# # Alpha = -1 : Neyman divergence-based

# # Beta = 2 : Euclid distance-based

# # Beta = 0 : Itakura-Saito divergence-based

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

# ####################################
# println("####### Simultaneous NMF (Command line) #######")

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

# ####################################
# println("####### Non-negative Tucker decomposition (Julia API) #######")
# # Alpha = 2 : Pearson divergence-based
# out_mu_alpha1 = ntd(input=joinpath(tmp, "X.zst"), alpha=2,
#   dim=3, numepoch=10, numbatch=10, logdir=tmp)

# @test size(out_mu_alpha1[1][1]) == (30, 3)
# @test size(out_mu_alpha1[1][2]) == (40, 3)
# @test size(out_mu_alpha1[1][3]) == (50, 3)
# @test size(out_mu_alpha1[2][1]) == (3, )
# @test size(out_mu_alpha1[2][2]) == (3, )
# @test size(out_mu_alpha1[2][3]) == (3, )
# @test size(out_mu_alpha1[3][1]) == (300, 3)
# @test size(out_mu_alpha1[3][2]) == (300, 3)
# @test size(out_mu_alpha1[3][3]) == (300, 3)
# @test size(out_mu_alpha1[4][1]) == ()
# @test size(out_mu_alpha1[4][2]) == ()
# @test size(out_mu_alpha1[4][3]) == ()
# @test size(out_mu_alpha1[5][1]) == ()
# @test size(out_mu_alpha1[5][2]) == ()
# @test size(out_mu_alpha1[5][3]) == ()
# @test size(out_mu_alpha1[6][1]) == ()
# @test size(out_mu_alpha1[6][2]) == ()
# @test size(out_mu_alpha1[6][3]) == ()
# @test size(out_mu_alpha1[7][1]) == ()
# @test size(out_mu_alpha1[7][2]) == ()
# @test size(out_mu_alpha1[7][3]) == ()
# @test size(out_mu_alpha1[8]) == ()

# # Alpha = 1 : KL divergence-based

# # Alpha = 0.5 : Hellinger divergence-based

# # Alpha = -1 : Neyman divergence-based

# # Beta = 2 : Euclid distance-based

# # Beta = 0 : Itakura-Saito divergence-based

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

# ####################################
# println("####### Non-negative Tucker decomposition (Command line) #######")

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


text:
using Rsvg
using Cairo
using Gtk

using Base.Test

include("test.jl")

pkg_dir = Pkg.dir("Rsvg")

@printf("\nTest: dimension of known images")

d = test_get_dimension(joinpath(pkg_dir,"data","mulberry.svg"));
@test (d.height == 512) & (d.width == 513)

d = test_get_dimension(joinpath(pkg_dir,"data","diag.svg"));
@test (d.height == 400) & (d.width == 400)

d = test_get_dimension(joinpath(pkg_dir,"data","lotus.svg"));
@test (d.height == 720) & (d.width == 576)

d = test_get_dimension(joinpath(pkg_dir,"data","star.svg"));
@test (d.height == 198) & (d.width == 224)

@printf("\nTest: render to png")

f = tempname() * ".png"
test_render_to_png(joinpath(pkg_dir,"data","mulberry.svg"),f);
@test stat(f).size > 0

f = tempname() * ".png"
test_render_to_png(joinpath(pkg_dir,"data","diag.svg"),f);
@test stat(f).size > 0

f = tempname() * ".png"
test_render_to_png(joinpath(pkg_dir,"data","lotus.svg"),f);
@test stat(f).size > 0

f = tempname() * ".png"
test_render_to_png(joinpath(pkg_dir,"data","star.svg"),f);
@test stat(f).size > 0

@printf("\nTest: render string to png")

f = tempname() * ".png"
test_render_string_to_png(f);
@test stat(f).size > 0

@printf("\nTest: roundtrip, render svg to svg")

f = tempname() * ".svg"
test_roundtrip(joinpath(pkg_dir,"data","lotus.svg"),f);
@test stat(f).size > 0
d = test_get_dimension(f);
@test (d.height == 720) & (d.width == 576)

@printf("\n");


text:
using BinaryBuilder

name = "Glib"
version = v"2.68.3"

# Collection of sources required to build Glib
sources = [
    ArchiveSource("https://ftp.gnome.org/pub/gnome/sources/glib/$(version.major).$(version.minor)/glib-$(version).tar.xz",
                  "e7e1a3c20c026109c45c9ec4a31d8dcebc22e86c69486993e565817d64be3138"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/glib-*/

if [[ "${target}" == *-freebsd* ]]; then
    # Our FreeBSD libc has `environ` as undefined symbol, so the linker will
    # complain if this symbol is used in the built library, even if this won't
    # be a problem at runtim.  This flag allows having undefined symbols.
    MESON_FLAGS=(-Db_lundef=false)

    # Adapt patch relative to `xattr` from
    # http://cvsweb.netbsd.org/bsdweb.cgi/pkgsrc/devel/glib2/patches/patch-meson.build?rev=1.2&content-type=text/x-cvsweb-markup.
    # Quoting the comment:
    #     Don't fail if getxattr is not available. The code is already ready
    #     for this case with some small configure changes.
    atomic_patch -p1 ../patches/freebsd-have_xattr.patch
fi

mkdir build_glib && cd build_glib
meson --cross-file="${MESON_TARGET_TOOLCHAIN}" \
    --buildtype=release \
    -Dman=false \
    -Diconv=external \
    "${MESON_FLAGS[@]}" \
    ..
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgio-2", "libgio-2.0"], :libgio),
    LibraryProduct(["libglib-2", "libglib-2.0"], :libglib),
    LibraryProduct(["libgmodule-2", "libgmodule-2.0"], :libgmodule),
    LibraryProduct(["libgobject-2", "libgobject-2.0"], :libgobject),
    LibraryProduct(["libgthread-2", "libgthread-2.0"], :libgthread),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Host gettext needed for "msgfmt"
    HostBuildDependency("Gettext_jll"),
    Dependency("Libiconv_jll"),
    Dependency("Libffi_jll", v"3.2.2"; compat="~3.2.2"),
    # Gettext is only needed on macOS, as far as I could see
    Dependency("Gettext_jll", v"0.21.0"; compat="=0.21.0"),
    Dependency("PCRE_jll"),
    Dependency("Zlib_jll"),
    Dependency("Libmount_jll"),
]

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

max_stars_repo_name: cometscome/Yggdrasil

max_issues_repo_name: cometscome/Yggdrasil

max_forks_repo_name: cometscome/Yggdrasil


text:
#-------------------------------------------------------------------------------
# Ply file IO functionality

#--------------------------------------------------
# Header IO

@enum Format Format_ascii Format_binary_little Format_binary_big

function ply_type(type_name)
    if     type_name == "char"   || type_name == "int8";    return Int8
    elseif type_name == "short"  || type_name == "int16";   return Int16
    elseif type_name == "int"    || type_name == "int32" ;  return Int32
    elseif type_name == "int64";                            return Int64
    elseif type_name == "uchar"  || type_name == "uint8";   return UInt8
    elseif type_name == "ushort" || type_name == "uint16";  return UInt16
    elseif type_name == "uint"   || type_name == "uint32";  return UInt32
    elseif type_name == "uint64";                           return UInt64
    elseif type_name == "float"  || type_name == "float32"; return Float32
    elseif type_name == "double" || type_name == "float64"; return Float64
    else
        error("type_name $type_name unrecognized/unimplemented")
    end
end

ply_type_name(::Type{UInt8})    = "uint8"
ply_type_name(::Type{UInt16})   = "uint16"
ply_type_name(::Type{UInt32})   = "uint32"
ply_type_name(::Type{UInt64})   = "uint64"
ply_type_name(::Type{Int8})     = "int8"
ply_type_name(::Type{Int16})    = "int16"
ply_type_name(::Type{Int32})    = "int32"
ply_type_name(::Type{Int64})    = "int64"
ply_type_name(::Type{Float32})  = "float32"
ply_type_name(::Type{Float64})  = "float64"

ply_type_name(::UInt8)    = "uint8"
ply_type_name(::UInt16)   = "uint16"
ply_type_name(::UInt32)   = "uint32"
ply_type_name(::UInt64)   = "uint64"
ply_type_name(::Int8)     = "int8"
ply_type_name(::Int16)    = "int16"
ply_type_name(::Int32)    = "int32"
ply_type_name(::Int64)    = "int64"
ply_type_name(::Float32)  = "float32"
ply_type_name(::Float64)  = "float64"

const PlyNativeType = Union{UInt8,UInt16,UInt32,UInt64,Int8,Int16,Int32,Int64,Float32,Float64}

ply_type_name(A::AbstractArray{T}) where {T<:PlyNativeType} = ply_type_name(T)

ply_type_name(A::AbstractArray) = !isempty(A) ? ply_type_name(A[1]) :
                                  error("Unknown ply element type name for empty array of type $(typeof(A))")

const _host_is_little_endian = (ENDIAN_BOM == 0x04030201)

function read_header(ply_file)
    firstline = readline(ply_file)
    if firstline != "ply"
        throw(ErrorException("Expected \"ply\" header, got \"$firstline\""))
    end
    element_name = ""
    element_numel = 0
    element_props = Vector{AbstractVector}()
    elements = PlyElement[]
    comments = PlyComment[]
    format = nothing
    while true
        line = readline(ply_file)
        if line == "" && eof(ply_file)
            throw(ErrorException("Unexpected end of file reading ply header"))
        elseif line == "end_header"
            break
        elseif startswith(line, "comment")
            push!(comments, PlyComment(strip(line[8:end]), false, length(elements)+1))
        elseif startswith(line, "obj_info")
            push!(comments, PlyComment(strip(line[9:end]), true, length(elements)+1))
        else
            tokens = split(line)
            length(tokens) > 2 || throw(ErrorException("Bad ply header, line: \"$line\""))
            if tokens[1] == "format"
                length(tokens) == 3 || throw(ErrorException("Bad ply header, line: \"$line\""))
                _, format_type, format_version = tokens
                if format_version != "1.0"
                    throw(ErrorException("Expected ply version 1.0, got $format_version"))
                end
                format = format_type == "ascii"                ? Format_ascii :
                        format_type == "binary_little_endian" ? Format_binary_little :
                        format_type == "binary_big_endian"    ? Format_binary_big :
                        error("Unknown ply format $format_type")
            elseif tokens[1] == "element"
                if !isempty(element_name)
                    push!(elements, PlyElement(element_name, element_numel, element_props))
                    element_props = Vector{AbstractVector}()
                end
                length(tokens) == 3 || throw(ErrorException("Bad ply header, line: \"$line\""))
                _, element_name, element_numel = tokens
                element_numel = parse(Int,element_numel)
            elseif tokens[1] == "property"
                !isempty(element_name) || throw(ErrorException("Bad ply header: property before first element"))
                if tokens[2] == "list"
                    length(tokens) == 5 || throw(ErrorException("Bad ply header, line: \"$line\""))
                    count_type_name, type_name, prop_name = tokens[3:end]
                    count_type = ply_type(count_type_name)
                    type_ = ply_type(type_name)
                    push!(element_props, ListProperty(prop_name, ply_type(count_type_name), ply_type(type_name)))
                else
                    length(tokens) == 3 || throw(ErrorException("Bad ply header, line: \"$line\""))
                    type_name, prop_name = tokens[2:end]
                    push!(element_props, ArrayProperty(prop_name, ply_type(type_name)))
                end
            else
                throw(ErrorException("Bad ply header, line: \"$line\""))
            end
        end
    end
    if !isempty(element_name)
        push!(elements, PlyElement(element_name, element_numel, element_props))
    end
    elements, format, comments
end

function write_header_field(stream::IO, prop::ArrayProperty)
    println(stream, "property $(ply_type_name(prop.data)) $(prop.name)")
end

function write_header_field(stream::IO, prop::ArrayProperty{T,Names}) where {T,Names<:PropNameList}
    for n in prop.name
        println(stream, "property $(ply_type_name(prop.data)) $(n)")
    end
end

function write_header_field(stream::IO, prop::ListProperty{S}) where {S}
    println(stream, "property list $(ply_type_name(S)) $(ply_type_name(prop.data)) $(prop.name)")
end

function write_header_field(stream::IO, comment::PlyComment)
    prefix = comment.obj_info ? "obj_info " : "comment "
    println(stream, prefix, comment.comment)
end

function write_header(ply, stream::IO, ascii)
    println(stream, "ply")
    if ascii
        println(stream, "format ascii 1.0")
    else
        endianness = _host_is_little_endian ? "little" : "big"
        println(stream, "format binary_$(endianness)_endian 1.0")
    end
    commentidx = 1
    for (elemidx,element) in enumerate(ply.elements)
        while commentidx <= length(ply.comments) && ply.comments[commentidx].location == elemidx
            write_header_field(stream, ply.comments[commentidx])
            commentidx += 1
        end
        println(stream, "element $(element.name) $(length(element))")
        for property in element.properties
            write_header_field(stream, property)
        end
    end
    while commentidx <= length(ply.comments)
        write_header_field(stream, ply.comments[commentidx])
        commentidx += 1
    end
    println(stream, "end_header")
end

#-------------------------------------------------------------------------------
# ASCII IO for properties and elements

function parse_ascii(::Type{T}, io::IO) where {T}
    # FIXME: sadly unbuffered, will probably have terrible performance.
    buf = UInt8[]
    while !eof(io)
        c = read(io, UInt8)
        if c == UInt8(' ') || c == UInt8('\t') || c == UInt8('\r') || c == UInt8('\n')
            if !isempty(buf)
                break
            end
        else
            push!(buf, c)
        end
    end
    parse(T, String(buf))
end

function read_ascii_value!(stream::IO, prop::ArrayProperty{T}, index) where {T}
    prop.data[index] = parse_ascii(T, stream)
end
function read_ascii_value!(stream::IO, prop::ListProperty{S,T}, index) where {S,T}
    N = parse_ascii(S, stream)
    prop.start_inds[index+1] = prop.start_inds[index] + N
    for i=1:N
        push!(prop.data, parse_ascii(T, stream))
    end
end

#-------------------------------------------------------------------------------
# Binary IO for properties and elements

#--------------------------------------------------
# property IO
function read_binary_value!(stream::IO, prop::ArrayProperty{T}, index) where {T}
    prop.data[index] = read(stream, T)
end
function read_binary_value!(stream::IO, prop::ListProperty{S,T}, index) where {S,T}
    N = read(stream, S)
    prop.start_inds[index+1] = prop.start_inds[index] + N
    inds = read!(stream, Vector{T}(undef, Int(N)))
    append!(prop.data, inds)
end

function write_binary_value(stream::IO, prop::ArrayProperty, index)
    write(stream, prop.data[index])
end
function write_binary_value(stream::IO, prop::ListProperty{S}, index) where {S}
    len = prop.start_inds[index+1] - prop.start_inds[index]
    write(stream, convert(S, len))
    esize = sizeof(eltype(prop.data))
    unsafe_write(stream, pointer(prop.data) + esize*(prop.start_inds[index]-1), esize*len)
end

function write_ascii_value(stream::IO, prop::ListProperty, index)
    print(stream, prop.start_inds[index+1] - prop.start_inds[index], ' ')
    for i = prop.start_inds[index]:prop.start_inds[index+1]-1
        if i != prop.start_inds[index]
            write(stream, ' ')
        end
        print(stream, prop.data[i])
    end
end
function write_ascii_value(stream::IO, prop::ArrayProperty, index)
    print(stream, prop.data[index])
end
function write_ascii_value(stream::IO, prop::ArrayProperty{<:AbstractArray}, index)
    p = prop.data[index]
    for i = 1:length(p)
        if i != 1
            write(stream, '\t')
        end
        print(stream, p[i])
    end
end

#--------------------------------------------------
# Batched element IO

# Read/write values for an element as binary.  We codegen a version for each
# number of properties so we can unroll the inner loop to get type inference
# for individual properties.  (Could this be done efficiently by mapping over a
# tuple of properties?  Alternatively a generated function would be ok...)
for numprop=1:16
    propnames = [Symbol("p$i") for i=1:numprop]
    @eval function write_binary_values(stream::IO, elen, $(propnames...))
        for i=1:elen
            $([:(write_binary_value(stream, $(propnames[j]), i)) for j=1:numprop]...)
        end
    end
    @eval function read_binary_values!(stream::IO, elen, $(propnames...))
        for i=1:elen
            $([:(read_binary_value!(stream, $(propnames[j]), i)) for j=1:numprop]...)
        end
    end
end
# Fallback for large numbers of properties
function write_binary_values(stream::IO, elen, props...)
    for i=1:elen
        for p in props
            write_binary_value(stream, property, i)
        end
    end
end
function read_binary_values!(stream::IO, elen, props...)
    for i=1:elen
        for p in props
            read_binary_value!(stream, p, i)
        end
    end
end

# Optimization: special cases for a single array property within an element
function write_binary_values(stream::IO, elen, prop::ArrayProperty)
    write(stream, prop.data)
end
function read_binary_values!(stream::IO, elen, prop::ArrayProperty)
    read!(stream, prop.data)
end

# Optimization: For properties with homogeneous type, shuffle into a buffer
# matrix before a batch-wise call to write().  This is a big speed improvement
# for elements constructed of simple arrays with homogenous type -
# serialization speed generally seems to be limited by the many individual
# calls to write() with small buffers.
function write_binary_values(stream::IO, elen, props::ArrayProperty{T}...) where {T}
    batchsize = 100
    numprops = length(props)
    buf = Matrix{T}(numprops, batchsize)
    for i=1:batchsize:elen
        thisbatchsize = min(batchsize, elen-i+1)
        for j=1:numprops
            buf[j,1:thisbatchsize] = props[j].data[i:i+thisbatchsize-1]
        end
        unsafe_write(stream, pointer(buf), sizeof(T)*numprops*thisbatchsize)
    end
end

#-------------------------------------------------------------------------------
# High level IO for complete files

"""
    load_ply(file)

Load data from a ply file and return a `Ply` datastructure.  `file` may either
be a file name or an open stream.
"""
function load_ply(io::IO)
    elements, format, comments = read_header(io)
    if format != Format_ascii
        if _host_is_little_endian && format != Format_binary_little
            error("Reading big endian ply on little endian host is not implemented")
        elseif !_host_is_little_endian && format != Format_binary_big
            error("Reading little endian ply on big endian host is not implemented")
        end
    end
    for element in elements
        for prop in element.properties
            resize!(prop, length(element))
        end
        if format == Format_ascii
            for i = 1:length(element)
                for prop in element.properties
                    read_ascii_value!(io, prop, i)
                end
            end
        else # format == Format_binary_little
            read_binary_values!(io, length(element), element.properties...)
        end
    end
    Ply(elements, comments)
end

function load_ply(file_name::AbstractString)
    open(file_name, "r") do fid
        load_ply(fid)
    end
end

"""
    save_ply(ply::Ply, file; [ascii=false])

Save data from `Ply` data structure into `file` which may be a filename or an
open stream.  The file will be native endian binary, unless the keyword
argument `ascii` is set to `true`.
"""
function save_ply(ply, stream::IO; ascii::Bool=false)
    write_header(ply, stream, ascii)
    for element in ply
        if ascii
            for i=1:length(element)
                for (j,property) in enumerate(element.properties)
                    if j != 1
                        write(stream, '\t')
                    end
                    write_ascii_value(stream, property, i)
                end
                println(stream)
            end
        else # binary
            write_binary_values(stream, length(element), element.properties...)
        end
    end
end

function save_ply(ply, file_name::AbstractString; kwargs...)
    open(file_name, "w") do fid
        save_ply(ply, fid; kwargs...)
    end
end


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Combine a feature along with the number of scenario successes and failures.
"""
struct FeatureSuccessAndFailure
    feature::Feature
    n_success::UInt
    n_failure::UInt
end

"""
Accumulate results from executed features as they are being executed. Keep track of whether the
total run is a success of a failure.
"""
mutable struct ResultAccumulator
    isaccumsuccess::Bool
    features::Vector{FeatureSuccessAndFailure}
    errors::Vector{Tuple{String, Gherkin.BadParseResult{Feature}}}

    ResultAccumulator() = new(true, [], [])
end

"""
    accumulateresult!(acc::ResultAccumulator, result::FeatureResult)

Check for success or failure in this feature result and update the accumulator accordingly.
"""
function accumulateresult!(acc::ResultAccumulator, result::FeatureResult)
    n_success::UInt = 0
    n_failure::UInt = 0
    # Count the number of successes and failures for each step in each scenario. A scenario is
    # successful if all its steps are successful.
    for scenarioresult in result.scenarioresults
        arestepssuccessful = [issuccess(step) for step in scenarioresult.steps]
        arebackgroundstepsgood = [issuccess(result) for result in scenarioresult.backgroundresult]
        isscenariosuccessful = all(arestepssuccessful) && all(arebackgroundstepsgood)
        if isscenariosuccessful
            n_success += 1
        else
            n_failure += 1
        end

        acc.isaccumsuccess &= isscenariosuccessful
    end

    push!(acc.features, FeatureSuccessAndFailure(result.feature, n_success, n_failure))
end

"""
    accumulateresult!(acc::ResultAccumulator, parsefailure::Gherkin.BadParseResult{Feature})

A feature file could not be parsed properly. Record the error.
"""
const BadParseType = Union{Gherkin.BadParseResult{Feature}, Gherkin.Experimental.BadParseResult{Feature}}
function accumulateresult!(acc::ResultAccumulator, parsefailure::BadParseType, featurefile::String)
    push!(acc.errors, (featurefile, parsefailure))
    acc.isaccumsuccess = false
end

"""
    issuccess(acc::ResultAccumulator)

True if all scenarios in all accumulated features are successful.
"""
issuccess(acc::ResultAccumulator) = acc.isaccumsuccess

"""
    featureresults(accumulator::ResultAccumulator)

A public getter for the results of all features.
"""
function featureresults(accumulator::ResultAccumulator)
    accumulator.features
end

"""
    isempty(r::ResultAccumulator) :: Bool

True if no results have been accumulated, false otherwise.
"""
Base.isempty(r::ResultAccumulator) :: Bool = isempty(r.features)


text:
mutable struct TemporalRegularizationOperator{V<:AbstractMatrix}
  S::V
  Φ::V
  N::Int
  M::Int
  L::Int # Number of all raw data frames
  idxFG::Vector{Int} # Sampling indices of foreground (sorted)
  idxBG::Vector{Int} # Sampling indices of background (sorted)
  idxCoeffsFG::Vector{Vector{Int}}
  idxCoeffsBG::Vector{Vector{Int}}
  coeffsFG::Vector{Vector{Float32}}
  coeffsBG::Vector{Vector{Float32}}
end

function reconstructionTempReg(S, bSF::Union{T,Vector{T}}, bMeas::MPIFile, freq::Array, grid;
  frames = nothing, bEmpty = nothing, emptyMeas= bEmpty, bgFrames = 1, numAverages = 1,
  bgDict = nothing, spectralLeakageCorrection=false, λ = 0.001, β = 0.01, solver = "kaczmarz",
  bgCorrectionInternal=false, interpMeth=:linear, idxFG=nothing, idxBG=nothing, kargs...) where {T<:MPIFile}

  bgCorrection = emptyMeas != nothing ? true : false

  @debug "Loading emptymeas ..."
  if emptyMeas != nothing
      uEmpty = getMeasurementsFD(emptyMeas, frequencies=freq, frames=bgFrames, numAverages=length(bgFrames),
       spectralLeakageCorrection=spectralLeakageCorrection, bgCorrection=bgCorrectionInternal)
  end

  frames == nothing && (frames = 1:acqNumFrames(bMeas))

  u = getMeasurementsFD(bMeas, frequencies=freq, frames=frames, numAverages=numAverages, spectralLeakageCorrection=spectralLeakageCorrection,
                          bgCorrection=bgCorrectionInternal)

  L = -fld(-length(frames),numAverages) # number of tomograms to be reconstructed

  idxFG == nothing && (idxFG = collect(1:L))
  idxBG == nothing && (idxBG = collect(1:L))

  if emptyMeas != nothing
    u = u .- uEmpty
  end

  @info size(u)  size(S) size(bgDict) L

  #initialize output
  #image = initImage(bSF,bMeas,L,numAverages,grid,false)

  if sum(abs.(λ)) > 0
    trace = calculateTraceOfNormalMatrix(S,nothing)
    λ *= trace / size(S,1)
  end

  Op = TemporalRegularizationOperator(Float32(1/(sqrt(λ)))*S, Float32(1/(sqrt(β)))*bgDict, L;
                                      interpMeth, idxFG, idxBG )

  global Top = Op

  c = reconstructionTempReg(Op, u; λ, β, solver, shape=shape(grid), kargs...)

  cInterp = similar(c, size(c,1), Op.L)
  for l=1:Op.L
    for κ=1:length(Op.idxCoeffsFG[l])
      cInterp[:,l] .+= Op.coeffsFG[l][κ] * c[:, Op.idxCoeffsFG[l][κ]]
    end
  end

  image = initImage(bSF, bMeas, Op.L, numAverages, grid, false)

  @info size(image)

  image[:] .= vec(cInterp)

  @info size(image)

  return image
end

global Top = nothing

"""
Joint reconstruction of MPI signal and background. Implemented in a low-level
fashion
"""
function reconstructionTempReg(Op::TemporalRegularizationOperator, u::Array;
                        β = 0.1, λ=0.001, solver = "kaczmarz", enforceReal=false, enforcePositive=false,
                        relativeLambda=true, backgroundCoefficients = nothing, kargs...)

  MSub, NSub, Q, J, Γ = getSizes(Op)

  u = reshape(u, :, Op.L)
  c = zeros(NSub, J)

  constraintMask = zeros(Bool, Op.N)
  constraintMask[1:(NSub*J)] .= 1

  solv = createLinearSolver(solver, Op; λ=1.0, constraintMask=constraintMask,
                            enforceReal=enforceReal,
			                      enforcePositive=enforcePositive, kargs...)

  @info "SIZES" size(Op) size(u) size(reshape(u,:,1))

  y = solve(solv, vec(u))
  c = real.( reshape(y[1:(NSub*J),:],NSub,J) ./ sqrt(λ) )

  #if backgroundCoefficients != nothing
  #  append!(backgroundCoefficients, vec(y[(N+1):end,:] ./ sqrt(β)))
  #end

  return c
end

getSizes(Op::TemporalRegularizationOperator) =
      Op.M÷Op.L, size(Op.S,2), size(Op.Φ,2), length(Op.idxFG), length(Op.idxBG)

Base.eltype(TempOp::TemporalRegularizationOperator) = eltype(TempOp.S)

function size(TempOp::TemporalRegularizationOperator, i::Int)
  if i==2
    return TempOp.N
  elseif i==1
    return TempOp.M
  else
    error("bounds error")
  end
end
size(TempOp::TemporalRegularizationOperator) = (TempOp.M,TempOp.N)
length(TempOp::TemporalRegularizationOperator) = TempOp.M*TempOp.N

function TemporalRegularizationOperator(S, Φ, L; interpMeth=:linear, idxFG=1:L, idxBG=1:L )

  # check if first and last frame are sampled
  if idxFG[1] != 1 || idxFG[end] != L ||
    idxBG[1] != 1 || idxBG[end] != L
    error("idxFG or idxBG do not contain first and last frame!")
  end

  idxCoeffsFG, idxCoeffsBG, coeffsFG, coeffsBG =
       calculateInterpolatingCoefficients(L, idxFG, idxBG, interpMeth)

  N = size(S,2)*length(idxFG)+size(S,2)*length(idxBG)
  M = size(S,1)*L

  @info "gluck" N M size(S)

  return TemporalRegularizationOperator(S, Φ, N, M, L, idxFG, idxBG,
                           idxCoeffsFG, idxCoeffsBG, coeffsFG, coeffsBG)
end

function calculateInterpolatingCoefficients(L, idxFG, idxBG, interpMeth)
  idxCoeffsFG = Vector{Vector{Int}}(undef,0)
  idxCoeffsBG = Vector{Vector{Int}}(undef,0)
  coeffsFG = Vector{Vector{Float32}}(undef,0)
  coeffsBG = Vector{Vector{Float32}}(undef,0)

  if interpMeth==:linear
    function help(idx, idxCoeffs, coeffs)
      l = 1
      for j=1:length(idx)
        r = idx[j]:idx[min(j+1,end)]
        #push!(idxCoeffs, [idx[j]])
        push!(idxCoeffs, [j])
        push!(coeffs, [1.0])
        l += 1
        for q=2:length(r)-1
          α = (q-1)/(length(r)-1)
          #push!(idxCoeffs, [idx[j],idx[j+1]] )
          push!(idxCoeffs, [j,j+1] )
          push!(coeffs, [1-α,α] )
          l += 1
        end
      end
    end
    help(idxFG, idxCoeffsFG, coeffsFG)
    help(idxBG, idxCoeffsBG, coeffsBG)
  else
    error("Interpolation method $(interpMeth) is right now not supported")
  end
  return idxCoeffsFG, idxCoeffsBG, coeffsFG, coeffsBG
end

### The following is intended to use the standard kaczmarz method ###

function calculateTraceOfNormalMatrix(Op::TemporalRegularizationOperator, weights)
  return calculateTraceOfNormalMatrix(Op.S, weights)
end

setlambda(::TemporalRegularizationOperator, ::Any) = nothing

getSubIdxRows(k,M) = mod1(k,M), (k-1)÷M+1

function RegularizedLeastSquares.dot_with_matrix_row(Op::TemporalRegularizationOperator,
                                                     x::AbstractArray{T}, k::Integer) where T

  MSub, NSub, Q, J, Γ = getSizes(Op)

  #@info "Blub"
  #@info MSub, NSub, Q, J, Γ
  #@info size(x), k

  i,l = getSubIdxRows(k, MSub)
  xFG = reshape( view(x, 1:NSub*J), NSub, J)
  xBG = reshape( view(x, (NSub*J+1):(NSub*J+Q*Γ)), Q, Γ)

  #@info l
  #@info size(Op.idxCoeffsFG)
  #@info Op.idxCoeffsFG[l]

  tmpFG = Op.coeffsFG[l][1] * xFG[:, Op.idxCoeffsFG[l][1]]
  for κ=2:length(Op.idxCoeffsFG[l])
    tmpFG .+= Op.coeffsFG[l][κ] * xFG[:, Op.idxCoeffsFG[l][κ]]
  end

  tmpBG = Op.coeffsBG[l][1] * xBG[:, Op.idxCoeffsBG[l][1]]
  for κ=2:length(Op.idxCoeffsBG[l])
    tmpBG .+= Op.coeffsBG[l][κ] * xBG[:, Op.idxCoeffsBG[l][κ]]
  end

  tmp = zero(T)
  @simd  for n = 1:NSub
    @inbounds tmp += Op.S[i,n]*tmpFG[n]
  end
  @simd  for q = 1:Q
    @inbounds tmp += Op.Φ[i,q]*tmpBG[q]
  end
  return tmp
end

function RegularizedLeastSquares.kaczmarz_update!(Op::TemporalRegularizationOperator,
                                                  x::AbstractArray, k::Integer, β)
  MSub, NSub, Q, J, Γ = getSizes(Op)
  i, l = getSubIdxRows(k, MSub)
  xFG = reshape( view(x, 1:NSub*J), NSub, J)
  xBG = reshape( view(x, (NSub*J+1):(NSub*J+Q*Γ)), Q, Γ)

  for κ=1:length(Op.idxCoeffsFG[l])
    α = β * conj(Op.coeffsFG[l][κ])
    #@simd
    for n = 1:NSub
      #@inbounds
      xFG[n,Op.idxCoeffsFG[l][κ]] += α * conj(Op.S[i,n])
    end
  end

  for κ=1:length(Op.idxCoeffsBG[l])
    α = β * conj(Op.coeffsBG[l][κ])
    #@simd
    for q = 1:Q
     # @info q Q
      #@inbounds
      xBG[q,Op.idxCoeffsBG[l][κ]] += α * conj(Op.Φ[i,q])
    end
  end
  return
end

function RegularizedLeastSquares.initkaczmarz(Op::TemporalRegularizationOperator,λ,weights::Vector)
  T = typeof(real(Op.S[1]))
  denom = T[]
  rowindex = Int64[]
  MSub, NSub, Q, S, Γ = getSizes(Op)

  s² = [rownorm²(Op.S,i) for i=1:MSub]
  φ² = [rownorm²(Op.Φ,i) for i=1:MSub]

  for l=1:Op.L
    for i=1:MSub
      if s²[i] > 0 || φ²[i] > 0
        k = i+MSub*(l-1)
        push!(denom,1 / (s²[i]*sum(abs.(Op.coeffsFG[l]).^2) +
                         φ²[i]*sum(abs.(Op.coeffsBG[l]).^2) + λ))
        push!(rowindex,k)
      end
    end
  end

  denom, rowindex
end


text:
## data.jl Test some functionality of the Data type
## (c) 2015 David A. van Leeuwen
using Logging
using Distributed

@info("Testing Data")

using JLD
for i = 1:10
    save("$i.jld", "data", randn(10000,3))
end
x = Matrix{Float64}[load("$i.jld", "data") for i=1:10]

using GaussianMixtures

g = rand(GMM, 2, 3)
d = Data(x)
dd = Data(["$i.jld" for i=1:10], Float64)

f1(gmm, data) = GaussianMixtures.dmapreduce(x->stats(gmm, x), +, data)
f2(gmm, data) = reduce(+, map(x->stats(gmm, x), data))

sleep(1)
println(f2(g,dd))

import Base.isapprox
isapprox(a::Array, b::Array) = size(a)==size(b) && all(i->isapprox(a[i], b[i]), 1:length(a))
isapprox(a::Tuple, b::Tuple) = all(Bool[isapprox(x,y) for (x,y) in zip(a,b)])

s = stats(g, collect(d))

@assert isapprox(s, f1(g,d))
@assert isapprox(s, f1(g,dd))
@assert isapprox(s, f2(g,d))
@assert isapprox(s, f2(g,dd))

@assert isapprox(s, stats(g,d))
@assert isapprox(s, stats(g,dd))
@assert isapprox(s, stats(g,d, parallel=true))
@assert isapprox(s, stats(g,dd, parallel=true))

rmprocs(workers())


text:
################################################################################
######################### Wind Calculations ####################################
################################################################################
#**
function wndF_wndPrf(profiles)
    wp=Array{Float32,1}()
    dummy=Array{Float32,1}()
    dummy=profiles[1]
    for i=2:length(profiles)
        dummy=dummy+profiles[i]
    end
    wp=dummy./length(profiles)

    wnd=wind()
    mx=findmax(wp)[1]
    ord=reverse(sort(wp))./mx
    wndF_conEng(ord,wnd)#constraint energy calc
    wndF_ldLss(ord, wnd)#calculates loss factor
    wnd
    return wnd
end
#loss factor
#**
function wndF_ldLss(div, wind)
  #wind.lf=(sum(div))*0.85/length(div)#saves loss factor, 0.85 is wake penalization
  wind.lf=(sum(div))/length(div)#saves loss factor, 0.85 is wake penalization
  #loss factor/llf formula ref: Guidelines on the calculation and use of loss factors Te Mana Hiko Electricity Authority
  #0.85 for wake effect from Evaluation of the wind direction uncertainty and its impact on wake modeling at the Horns Rev offshore wind farm
  #M. Gaumond  P.‐E. Réthoré  S. Ott  A. Peña  A. Bechmann  K. S. Hansen
  llf=0.0
  for pu in div
    #llf=llf+(pu*0.85)^2
    llf=llf+pu^2
  end
  wind.delta=llf/length(div)#saves load loss factor
end

#Finds Constraind energy
#**
function  wndF_conEng(graph,wnd)
    #create sized arrays
    ce=Float32[]
    for hr=1:length(graph)
        smPu=0
        for pu=1:hr
            smPu=smPu+graph[pu]
        end
        smPu=smPu-(graph[hr]*hr)
        push!(ce,smPu)
    end
    wnd.ce=deepcopy(ce)
    wnd.pu=deepcopy(graph)
    return nothing
end

#check
#**Totals wind profiles of set
function find_netWind(wndstring)
    if !(haskey(wind_module.wind_profs,wndstring))
        wind_sum=wind()
        wndstring=wndstring[1]=='_' ? wndstring[2:end] : wndstring;
        wnds=split(wndstring,"_")
        if (length(wnds)>0 && wnds[length(wnds)]!="")
            wind_sum.pu=zeros(Float32,8759)
            wind_sum.ce=zeros(Float32,8759)
            wind_sum.delta=0
            wind_sum.lf=0
            true_length=0
            for w in wnds
                if (length(wind_module.wind_profs[w].pu)>0)
                    wind_sum.pu=(wind_sum.pu.+wind_module.wind_profs[w].pu)
                    wind_sum.ce=(wind_sum.ce.+wind_module.wind_profs[w].ce)
                    wind_sum.delta=(wind_sum.delta+wind_module.wind_profs[w].delta)
                    wind_sum.lf=(wind_sum.lf+wind_module.wind_profs[w].lf)
                    true_length=true_length+1
                end
            end
            if true_length!=0
                wind_sum.pu=wind_sum.pu./true_length
                wind_sum.ce=wind_sum.ce./true_length
                wind_sum.delta=wind_sum.delta/true_length
                wind_sum.lf=wind_sum.lf/true_length
            end
        end
        #println("storing: "*wndstring)
        wind_module.save_wind4_module(wind_sum,wndstring)
    else
    end
    return wndstring
end


text:
function TwoPhase(Kx::Array{T,3}, ϕ::Array{T,3}, qinj::Tuple{T1, T1, T1}, qrate::Number, d::Tuple{T2, T2, T2}, time::Number, nt::Int; Ky=nothing, Kz=nothing, o=nothing) where {T, T1, T2}
    "Kx permeability, ϕ porosity, qinj injection coordinate [m], qrate injection rate [Mt/y]"

    if isnothing(Ky)
        Ky = Kx
    end
    if isnothing(Kz)
        Kz = Kx
    end
    if isnothing(o)
        o = zeros(Float32, 3)
    end

    WriteTxtFile(Kx; name="PERMX")
    WriteTxtFile(Ky; name="PERMY")
    WriteTxtFile(Kz; name="PERMZ")
    WriteTxtFile(ϕ; name="PORO")

    n = size(Kx)
    WriteDATAFile(qinj, qrate, n, d, o, time, nt;
    PERMXtxt="PERMX.txt", PERMYtxt="PERMY.txt", PERMZtxt="PERMZ.txt", POROtxt="PORO.txt", name="CO2SIMULATION")
    run(`flow CO2SIMULATION.DATA`)
    sat, p = ReadResults(; name="CO2SIMULATION")

    return sat, p
end

function TwoPhase(Kx::Array{T,2}, ϕ::Array{T,2}, qinj::Tuple{T1, T1}, qrate::Number, d::Tuple{T2, T2, T2}, time::Number, nt::Int; Ky=nothing, Kz=nothing, o=nothing) where {T, T1, T2}
    n = size(Kx)
    if isnothing(Ky)
        Ky = reshape(Kx,n[1],1,n[2])
    else
        Ky = reshape(Ky,n[1],1,n[2])
    end
    if isnothing(Kz)
        Kz = reshape(Kx,n[1],1,n[2])
    else
        Kz = reshape(Kz,n[1],1,n[2])
    end
    if isnothing(o)
        o = zeros(Float32, 3)
    else
        o = (o[1],0f0,o[2])
    end
    sat, p = TwoPhase(reshape(Kx,n[1],1,n[2]), reshape(ϕ,n[1],1,n[2]), (qinj[1],d[2]/2,qinj[2]), qrate, d, time, nt; Ky=Ky, Kz=Kz, o=o)
    sat = [dropdims(sat[i], dims=2) for i = 1:length(sat)]
    p = [dropdims(p[i], dims=2) for i = 1:length(p)]

    return sat, p
end

function ReadResults(;name="CO2SIMULATION")
    py"""
    from ecl.grid import EclGrid, EclRegion
    from ecl.eclfile import EclFile, EclRestartFile
    import matplotlib.pyplot as plt
    import numpy as np
    ###################################################################################################
    # Read eclipse files
    def readecl(name):
        # Load grid
        grid = EclGrid(name + ".EGRID")
        shape = (grid.nx, grid.ny, grid.nz)[::-1]
        # Load snapshots
        rst_file = EclRestartFile(grid, name + ".UNRST")
        sat = [rst_file['SGAS'][i].numpy_view().reshape(shape) for i in range(len(rst_file['SGAS']))]
        p = [rst_file['PRESSURE'][i].numpy_view().reshape(shape) for i in range(len(rst_file['PRESSURE']))]
        return sat, p
    """
    sat, p = py"readecl"(name)
    sat = permutedims.(sat, [[3,2,1] for i = 1:length(sat)])
    p = permutedims.(p, [[3,2,1] for i = 1:length(sat)])
    return sat, p
end

function ReadInitFile(;name="CO2SIMULATION")
    py"""
    from ecl.grid import EclGrid, EclRegion
    from ecl.eclfile import EclFile, EclRestartFile
    import matplotlib.pyplot as plt
    import numpy as np
    ###################################################################################################
    # Read eclipse files
    def readeclinit(name):
        # Load grid
        grid = EclGrid(name + ".EGRID")
        shape = (grid.nx, grid.ny, grid.nz)[::-1]
        init_file = EclFile(name + ".INIT")
        ϕ = init_file['PORO'][0].numpy_view().reshape(shape)
        K = init_file['PERMX'][0].numpy_view().reshape(shape)
        return K, ϕ
    """
    K, ϕ = py"readeclinit"(name)
    K = permutedims(K, [3,2,1])
    ϕ = permutedims(ϕ, [3,2,1])
    return K, ϕ
end

function WriteTxtFile(var::Array{T,3}; name="PERMX") where T

    println("Writing $(name).txt")
    f = open(name*".txt","w")

    write(f, name)
    write(f, "\n")
    global ct = 0
    for k = 1:size(var,3)
        for j = 1:size(var,2)
            for i = 1:size(var,1)
                global ct = ct + 1
                write(f, string(var[i,j,k]))
                write(f, " ")
                (ct%5 == 0) && write(f, "\n")
            end
        end
    end

    write(f, " /")
    close(f)

end

function WriteTxtFile(var::Array{T,2}; name="PERMX.txt") where T
    WriteTxtFile(reshape(var,size(var,1),1,size(var,2)); name=name)
end

function WriteDATAFile(qinj::Tuple{T1, T1, T1}, qrate::Number, n::Tuple{Int, Int, Int}, d::Tuple{T3, T3, T3}, o::Tuple{T2, T2, T2}, time::Number, nt::Int;
    PERMXtxt::String = "PERMX.txt", PERMYtxt::String = "PERMY.txt", PERMZtxt::String = "PERMZ.txt", POROtxt::String = "PORO.txt", name="CO2SIMULATION") where {T1, T2, T3}

    println("Writing $(name).txt")
    DIMENS = "DIMENS
$(n[1]) $(n[2]) $(n[3]) /"

    DX = "DX
    $(prod(n))*$(d[1]) /"

    DY = "DY
    $(prod(n))*$(d[2]) /"

    DZ = "DZ
    $(prod(n))*$(d[3]) /"

    TOPS = "TOPS
    $(prod(n))*$(o[3]) /"

    PERMX = "INCLUDE
    '$(PERMXtxt)'
/"

    PERMY = "INCLUDE
    '$(PERMYtxt)'
/"

    PERMZ = "INCLUDE
    '$(PERMZtxt)'
/"

    PORO = "INCLUDE
    '$(POROtxt)'
/"

    qgrid = Int.(round.((qinj.-o)./d))

    rate = qrate * 1f6 * 1f3 / 1.98 / 365

    WELL = "WELSPECS
    Injector I $(qgrid[1]) $(qgrid[2]) 0.0e+00 WATER 0 STD SHUT NO 0 SEG 0/
    /

    COMPDAT
    Injector $(qgrid[1]) $(qgrid[2]) $(qgrid[3]) $(qgrid[3]) OPEN -1 8.5107288246779274e+02 2.0e-01 -1.0 0 1* Y -1.0/
    /

    WCONINJE
    Injector GAS OPEN RATE $(rate) 3* 0 /
    /"

    TSTEP = "TSTEP
    $(nt)*$(time*365.25/nt)
    /
    "

    str =

"RUNSPEC

TITLE
    Simulation

$(DIMENS)

EQLDIMS
    1 100 50 1 50
/

TABDIMS
    1 1 40 20 2 20
/

WELLDIMS
    1 1 1 1
/

OIL
GAS
CO2STOR

METRIC
UNIFOUT
START
    1 'JAN' 2000
/
------------------------------------------------------
GRID

INIT

$(DX)

$(DY)

$(DZ)

$(TOPS)

$(PERMX)

$(PERMY)

$(PERMZ)

$(PORO)

------------------------------------------------------
------------------------------------------------------
PROPS

ROCK
    1.0e+01 1.6e-06
/

-- Brooks-Corey with coefficient 2.8 and
-- entry pressure 2.5kPa (Cavanagh. A. 2013)
SGOF
-- Column 1: gas saturation
-- Column 2: gas relative permeability
-- Column 3: oil relative permeability when oil, gas and connate water are present
-- Column 4: oil-gas capillary pressure
    0.0	0.0	1.0 	0.025
    0.1	0.0     0.740	0.026
    0.2	0.009	0.528	0.027
    0.3	0.030	0.359	0.028
    0.4	0.070	0.230	0.030
    0.5	0.136	0.136	0.032
    0.6	0.230	0.070	0.035
    0.7	0.359	0.030	0.038
    0.8	0.528	0.009	0.044
    0.9	0.740	0.000	0.057 /

SALINITY
    0.7 / --35-40g/l  -> 35-40g /kg -> 0.63-0.72 mol/g

---------------------------------------------------------
REGIONS
---------------------------------------------------------

------------------------------------------------------
SOLUTION

EQUIL
    8.6300e+02 1.465583e+02  5.050e+03 0.0 1.0e+02 0.0  1 0 0 /
/

RSVD
    $(o[end]) 0.0
    $(n[end]*d[end]+o[end]) 0.0
/

RTEMPVD
    800 32.7134
    1000 41
/
------------------------------------------------------
SUMMARY

-----------------------------------------------------
SCHEDULE
RPTSCHED
PRES SGAS RS WELLS
/
RPTRST
BASIC=1
/

$(WELL)

$(TSTEP)

END
"

    f = open(name*".DATA","w")

    write(f, str)
    close(f)
end

max_stars_repo_name: slimgroup/TwoPhaseSpeedTest

max_issues_repo_name: slimgroup/TwoPhaseSpeedTest

max_forks_repo_name: slimgroup/TwoPhaseSpeedTest


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@testset "PIBasis"  begin

##

using ACE1, Random
using Printf, Test, LinearAlgebra, JuLIP, JuLIP.Testing
using JuLIP: evaluate, evaluate_d
using ACE1.Testing: println_slim

##

@info("Basic test of PIBasis construction and evaluation")

ord = 5
maxdeg = 10
r0 = 1.0
rcut = 3.0
trans = PolyTransform(1, r0)
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
D = ACE1.SparsePSHDegree()
P1 = ACE1.BasicPSH1pBasis(Pr; species = :X, D = D)

dagbasis = ACE1.PIBasis(P1, ord, D, maxdeg)
basis = standardevaluator(dagbasis)

# check single-species
Nat = 15
Rs, Zs, z0 = ACE1.rand_nhd(Nat, Pr, :X)
AA = evaluate(basis, Rs, Zs, z0)
println_slim(@test(length(basis) == length(AA)))

# construct a basis with dag-evaluator, and check they are identical!!!
AAdag = evaluate(dagbasis, Rs, Zs, z0)
println_slim(@test(AA ≈ AAdag))

dAA = evaluate_d(basis, Rs, Zs, z0
   )
dAAdag = evaluate_d(dagbasis, Rs, Zs, z0)
println_slim(@test dAA ≈ dAAdag)

println_slim(@test all(JuLIP.Testing.test_fio(dagbasis)))

## check multi-species
maxdeg = 5
ord = 3
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
species = [:C, :O, :H]
P1 = ACE1.BasicPSH1pBasis(Pr; species = [:C, :O, :H], D = D)
dagbasis = ACE1.PIBasis(P1, ord, D, maxdeg)
basis = standardevaluator(dagbasis)
Rs, Zs, z0 = ACE1.rand_nhd(Nat, Pr, species)
AA = evaluate(basis, Rs, Zs, z0)
println(@test(length(basis) == length(AA)))
dAA = evaluate_d(basis, Rs, Zs, z0)
println(@test(length(basis) == length(AA)))

# construct a basis with dag-evaluator, and check they are identical!!!
AAdag = evaluate(dagbasis, Rs, Zs, z0)
println_slim(@test(AA ≈ AAdag))

dAA = evaluate_d(basis, Rs, Zs, z0)
dAAdag = evaluate_d(dagbasis, Rs, Zs, z0)
println_slim(@test dAA ≈ dAAdag)

println(@test all(JuLIP.Testing.test_fio(dagbasis)))
##

@info("check SparseDegree alias")
D_ = ACE1.SparseDegree()
println_slim(@test D == D_)

##

@info("Check several properties of PIBasis")
for species in (:X, :Si, [:C, :O, :H]), N = 1:5
   local AA, AAdag, dAA, dAAdag, dagbasis, basis, Rs, Zs, z0
   maxdeg = 7
   Nat = 15
   P1 = ACE1.BasicPSH1pBasis(Pr; species = species)
   dagbasis = ACE1.PIBasis(P1, N, D, maxdeg)
   basis = standardevaluator(dagbasis)
   @info("species = $species; N = $N; length = $(length(basis))")
   @info("test (de-)serialisation")
   # only require dagbasis to deserialize correctly
   println(@test all(JuLIP.Testing.test_fio(dagbasis)))
   @info("Check Permutation invariance")
   for ntest = 1:20
      Rs, Zs, z0 = ACE1.rand_nhd(Nat, Pr, species)
      p = randperm(length(Rs))
      print_tf(@test(evaluate(basis, Rs, Zs, z0) ≈
                     evaluate(basis, Rs[p], Zs[p], z0)))
   end
   println()
   @info("Check gradients")
   for ntest = 1:20
      Rs, Zs, z0 = ACE1.rand_nhd(Nat, Pr, species)
      AA = evaluate(basis, Rs, Zs, z0)
      dAA = evaluate_d(basis, Rs, Zs, z0)
      Us = [ rand(eltype(Rs)) .- 0.5 for _=1:length(Rs) ]
      dAA_dUs = transpose.(dAA) * Us
      errs = []
      for p = 2:12
         h = 0.1^p
         AA_h = evaluate(basis, Rs + h * Us, Zs, z0)
         dAA_h = (AA_h - AA) / h
         # @show norm(dAA_h - dAA_dUs, Inf)
         push!(errs, norm(dAA_h - dAA_dUs, Inf))
      end
      success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
      print_tf(@test success)
   end

   println()
   @info("Check Standard=DAG Evaluator")
   for ntest = 1:20
      Rs, Zs, z0 = ACE1.rand_nhd(Nat, Pr, species)
      AA = evaluate(basis, Rs, Zs, z0)
      AAdag = evaluate(dagbasis, Rs, Zs, z0)
      print_tf(@test AA ≈ AAdag)

      dAA = evaluate_d(basis, Rs, Zs, z0)
      dAAdag = evaluate_d(dagbasis, Rs, Zs, z0)
      print_tf(@test dAA ≈ dAAdag)
   end
   println()
end
println()

##

end


text:
struct Nil end

struct List{T}
    head::T
    tail::Union{List{T}, Nil}
end

# built a list from an array
mklist(array::AbstractArray{T}) where T =
    foldr(List{T}, array, init=Nil())

# implement the iteration protocol
Base.iterate(l::List) = iterate(l, l)
Base.iterate(::List, l::List) = l.head, l.tail
Base.iterate(::List, ::Nil) = nothing

max_stars_repo_name: ninjaaron/dispatching-design-patterns

max_issues_repo_name: ninjaaron/dispatching-design-patterns

max_forks_repo_name: ninjaaron/dispatching-design-patterns


text:
export AbstractAgent,
       update!, buffer, learner

import DataStructures:update!

"""
    AbstractAgent

Supertype of **agents**.
Usually, an agent needs to contain at least an [`AbstractLearner`](@ref) and an [`AbstractTurnBuffer`](@ref).

| Required Methods| Brief Description |
|:----------------|:------------------|
| `agent(obs)` | `agent`, an instance of an `AbstractAgent`, must be a functional object to receive an observation as input and return a pair of state and action (s => a) |
| `update!(agent, s, a, r, d, s′[, a′])` | Update the agent after an interaction with environment |
"""
abstract type AbstractAgent end

function update! end
buffertype(agent::AbstractAgent) = typeof(agent.buffer)

buffer(agent::AbstractAgent) = agent.buffer
learner(agent::AbstractAgent) = agent.learner