text:
"""
    C_iio_buffer_cancel(buffer)

Cancel all buffer operations.

# Parameters
- `buffer::Ptr{iio_buffer}` :The buffer for which operations should be canceled

This function cancels all outstanding buffer operations previously scheduled.
This means any pending `iio_buffer_push()` or `iio_buffer_refill()` operation will abort and return immediately,
any further invocations of these functions on the same buffer will return immediately with an error.

Usually `iio_buffer_push()` and `iio_buffer_refill()` will block until either all data has been transferred or a timeout occurs.
This can depending on the configuration take a significant amount of time.
`iio_buffer_cancel()` is useful to bypass these conditions if the buffer operation is supposed to be stopped in response to an external event (e.g. user input).

To be able to capture additional data after calling this function the buffer should be destroyed and then re-created.

This function can be called multiple times for the same buffer, but all but the first invocation will be without additional effect.

This function is thread-safe, but not signal-safe, i.e. it must not be called from a signal handler.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga0e42431688750313cfa077ab4f6e0282)
"""
function C_iio_buffer_cancel(buffer::Ptr{iio_buffer})
    ccall(
        (:iio_buffer_cancel, libIIO),
        Cvoid, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_destroy(buffer)

Destroy the given buffer.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# NOTE
After that function, the iio_buffer pointer shall be invalid.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#gaba58dc2780be63fead6f09397ce90d10)
"""
function C_iio_buffer_destroy(buffer::Ptr{iio_buffer})
    ccall(
        (:iio_buffer_destroy, libIIO),
        Cvoid, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_end(buffer)

Get the address that follows the last sample in a buffer.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# Returns
- A pointer corresponding to the address that follows the last sample present in the buffer

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#gab5300f917bbdfc5dafc093a60138f131)
"""
function C_iio_buffer_end(buffer::Ptr{iio_buffer})
    return ccall(
        (:iio_buffer_end, libIIO),
        Ptr{Cuchar}, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_first(buffer, channel)

Find the first sample of a channel in a buffer.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure
- `channel::Ptr{iio_channel}` : A pointer to an iio_channel structure

# Returns
- A pointer to the first sample found, or to the end of the buffer if no sample for the given channel is present in the buffer

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga000d2f4c8b72060db1c38ec905bf4156)
"""
function C_iio_buffer_first(buffer::Ptr{iio_buffer}, channel::Ptr{iio_channel})
    return ccall(
        (:iio_buffer_first, libIIO),
        Ptr{Cuchar}, (Ptr{iio_buffer}, Ptr{iio_channel}),
        buffer, channel
    );
end

"""
    C_iio_buffer_foreach_sample()

THIS IS A PLACEHOLDER. THE DOCUMENTATION BELOW IS ONLY A COPY/PASTE OF THE C DOCUMENTATION.

Call the supplied callback for each sample found in a buffer.

# Parameters
- buf : A pointer to an iio_buffer structure
- callback : A pointer to a function to call for each sample found
- data : A user-specified pointer that will be passed to the callback

# Returns
- number of bytes processed.

# NOTE
The callback receives four arguments:
- A pointer to the iio_channel structure corresponding to the sample,
- A pointer to the sample itself,
- The length of the sample in bytes,
- The user-specified pointer passed to `iio_buffer_foreach_sample`.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga810ec50155e82331b18ec71d3c507104)
"""
function C_iio_buffer_foreach_sample()
    return "PLACEHOLDER"
end

"""
    C_iio_buffer_get_data(buffer)

Retrieve a previously associated pointer of an iio_buffer structure.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# Returns
- The pointer previously associated if present, or NULL

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#gac110da795a50dc45fe998cace656329b)
"""
# can return null
function C_iio_buffer_get_data(buffer::Ptr{iio_buffer})
    return ccall(
        (:iio_buffer_get_data, libIIO),
        Ptr{Cuchar}, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_get_device(buffer)

Retrieve a pointer to the iio_device structure.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# Returns
- A pointer to an iio_device structure

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga42367567d47f501d1922d1b331cf64fb)
"""
function C_iio_buffer_get_device(buffer::Ptr{iio_buffer})
    return ccall(
        (:iio_buffer_get_device, libIIO),
        Ptr{iio_device}, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_get_poll_fd(buffer)

Get a pollable file descriptor.

Can be used to know when `iio_buffer_refill()` or `iio_buffer_push()` can be called

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# Returns
- On success, valid file descriptor
- On error, a negative errno code is returned

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga2ae96ee9f0748e55dfad996d6e9883f2)
"""
function C_iio_buffer_get_poll_fd(buffer::Ptr{iio_buffer})
    return ccall(
        (:iio_buffer_get_poll_fd, libIIO),
        Cint, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_push(buffer)

Send the samples to the hardware.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# Returns
- On success, the number of bytes written is returned
- On error, a negative errno code is returned

# NOTE
Only valid for output buffers

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#gae7033c625d128667a56cf482aa3149bd)
"""
function C_iio_buffer_push(buffer::Ptr{iio_buffer})
    return ccall(
        (:iio_buffer_push, libIIO),
        Cssize_t, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_push_partial(buffer, sample_count)

Send a given number of samples to the hardware.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure
- `samples_count::UInt` : The number of samples to submit

# Returns
- On success, the number of bytes written is returned
- On error, a negative errno code is returned

# NOTE
Only valid for output buffers

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga367b7368532aebb35a0d56bccc550570)
"""
function C_iio_buffer_push_partial(buffer::Ptr{iio_buffer}, sample_count::UInt)
    return ccall(
        (:iio_buffer_push_partial, libIIO),
        Cssize_t, (Ptr{iio_buffer}, Csize_t),
        buffer, sample_count
    );
end

"""
    C_iio_buffer_refill(buffer)

Fetch more samples from the hardware.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# Returns
- On success, the number of bytes read is returned
- On error, a negative errno code is returned

# NOTE
Only valid for input buffers

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#gac999e5244b5a2cbbca5ecaef8303a4ff)
"""
function C_iio_buffer_refill(buffer::Ptr{iio_buffer})
    return ccall(
        (:iio_buffer_refill, libIIO),
        Cssize_t, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_set_blocking_mode(buffer, blocking)

Make `iio_buffer_refill()` and `iio_buffer_push()` blocking or not.

After this function has been called with blocking == false, `iio_buffer_refill()` and `iio_buffer_push()` will return -EAGAIN if no data is ready.
A device is blocking by default.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure
- `blocking::Bool` : true if the buffer API should be blocking, else false

# Returns
- On success, 0
- On error, a negative errno code is returned

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#gadf834d825ece149886283bcb8c2a5466)
"""
function C_iio_buffer_set_blocking_mode(buffer::Ptr{iio_buffer}, blocking::Bool)
    return ccall(
        (:iio_buffer_set_blocking_mode, libIIO),
        Cint, (Ptr{iio_buffer}, Cuchar),
        buffer, blocking
    );
end

"""
    C_iio_buffer_set_data(buffer, data)

Associate a pointer to an iio_buffer structure.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure
- `data` : A pointer to the data to be associated. Must be castable to Ptr{Cuchar}.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga07f485e4e2de57c8c1cd0141611187dc)
"""
# Ptr{Cvoid} again
# also probably need to make sure data lifetime > buffer lifetime
function C_iio_buffer_set_data(buffer::Ptr{iio_buffer}, data)
    return ccall(
        (:iio_buffer_set_data, libIIO),
        Cvoid, (Ptr{iio_buffer}, Ptr{Cuchar}),
        buffer, data
    );
end

"""
    C_iio_buffer_start(buffer)

Get the start address of the buffer.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# Returns
- A pointer corresponding to the start address of the buffer

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga7fdacbfda79aa5120f34ea73ae2ea5ab)
"""
function C_iio_buffer_start(buffer::Ptr{iio_buffer})
    return ccall(
        (:iio_buffer_start, libIIO),
        Ptr{Cuchar}, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_buffer_step(buffer)

Get the step size between two samples of one channel.

# Parameters
- `buffer::Ptr{iio_buffer}` : A pointer to an iio_buffer structure

# Returns
- the difference between the addresses of two consecutive samples of one same channel

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#ga5532665a8776cec1c209d6cf8d0bb409)
"""
function C_iio_buffer_step(buffer::Ptr{iio_buffer})
    return ccall(
        (:iio_buffer_step, libIIO),
        Cssize_t, (Ptr{iio_buffer},),
        buffer
    );
end

"""
    C_iio_device_create_buffer(device, samples_count, cyclic)

Create an input or output buffer associated to the given device.

# Parameters
- `device::Ptr{iio_device}`: A pointer to an iio_device structure
- `samples_count::UInt` : The number of samples that the buffer should contain
- `cyclic::Bool` : If True, enable cyclic mode

# Returns
- On success, a pointer to an iio_buffer structure
- On error, if the assertions are enabled, throws an error.
- On error, if the assertions are disabled, returns NULL.

# NOTE
Channels that have to be written to / read from must be enabled before creating the buffer.

[libIIO documentation](https://analogdevicesinc.github.io/libiio/master/libiio/group__Buffer.html#gaea8067aca27b93a1260a0c563607a501)
"""
function C_iio_device_create_buffer(device::Ptr{iio_device}, samples_count::UInt, cyclic::Bool)
    @assert_null_pointer buffer = ccall(
        (:iio_device_create_buffer, libIIO),
        Ptr{iio_buffer}, (Ptr{iio_device}, Csize_t, Cuchar),
        device, samples_count, cyclic
    );
    return buffer;
end


text:
# This file defines staionary arm models

"""
    StationaryArmBase

`StationaryArmBase` is the concrete type for all the stationary arm models. Any
subtype should implement ```pull!()``` method
"""
abstract type StationaryArmBase <: ArmBase end

# Default method for interface
pull!( arm::StationaryArmBase )     = nothing
tick!( arm::StationaryArmBase )     = nothing
reset!( arm::StationaryArmBase )    = nothing
mean( arm::StationaryArmBase )      = Distributions.mean( arm.armDist )

"""
    Bernoulli( p )
`Bernoulli` creates a Bernoulli arm with success probability `p`.
"""
mutable struct Bernoulli <: StationaryArmBase
    armDist::Distributions.Bernoulli

    function Bernoulli( p::Real )
        new( Distributions.Bernoulli(p) )
    end
end

function pull!( arm::Bernoulli )
    return Distributions.rand(arm.armDist)
end

"""
    Beta( α, β )

`Beta` creates a Beta arm with parameters `α` and `β`.
"""
mutable struct Beta <: StationaryArmBase
    armDist::Distributions.Beta

    function Beta( α::Real, β::Real )
        new( Distributions.Beta( α, β ) )
    end
end

function pull!( arm::Beta )
    return Distributions.rand( arm.armDist )
end

"""
    Exponetial( θ )

`Expoenential` creates a Expoenentially distributed arm with mean reward `θ`
"""
mutable struct Exponential <: StationaryArmBase
    armDist::Distributions.Exponential

    function Exponential( λ::Real ) # λ is the mean parameter
        new( Distributions.Exponential(λ) )
    end
end

function pull!( arm::Exponential )
    return Distributions.rand( arm.armDist )
end

"""
    Normal( μ, σ )

Creates an arm with reward distributed as Normal random varaible with mean `μ`
and standard deviation `σ`.
"""
mutable struct Normal <: StationaryArmBase
    armDist::Distributions.Normal

    function Normal( μ::Real, σ::Real )
        new(
            Distributions.Normal( μ, σ )
        )
    end
end

function pull!( arm::Normal )
    return rand( arm.armDist )
end

ext: jl


text:
@testset "Kronecker sums" begin

    A = rand(4,4); IA = oneunit(A)
    B = rand(3,3); IB = oneunit(B)

    KS = A ⊕ B
    kronsum = kron(A, IB) + kron(IA, B)

    @test collect(KS) ≈ kronsum
    @test collect(KS) isa SparseMatrixCSC

    @test issquare(KS)

    C = rand(5,5); IC = oneunit(C)
    KS3 = A ⊕ B ⊕ C
    kronsum3 = kron(A,IB,IC) + kron(IA,B,IC) + kron(IA,IB,C)

    @test collect(KS3) ≈ kronsum3
    @test collect(KS3) isa SparseMatrixCSC

    @test order(KS) == 2
    @test order(KS3) == 3

    @test getmatrices(KS) == (A,B)

    @test getindex(KS,2,3) == kronsum[2,3]
    @test getindex(KS3,2,3) == kronsum3[2,3]

    D = rand(ComplexF64,6,6); ID = oneunit(D)

    @testset "Structure of sums" begin
        @test size(A ⊕ D) == size(A) .* size(D,1)
        @test size(B ⊕ C ⊕ D) == size(B) .* size(C) .* size(D)
        @test eltype(A ⊕ B) == Float64
        @test eltype(C ⊕ D) == ComplexF64
    end

    @testset "Basic linear algebra for sums" begin
        @test tr(KS) ≈ tr(kronsum)
        @test KS' == kronsum'
        @test transpose(KS) == transpose(kronsum)
        @test conj(KS) == conj(kronsum)
    end

    A = rand(10,10); B = rand(10,10); V = Diagonal(rand(10))
    @testset "Vec trick for sums" begin
        @test (A ⊕ B) * vec(V) == vec(B*V + V*transpose(A))
    end

    A = rand(3, 3); IA = oneunit(A)
    B = rand(4, 4); IB = oneunit(B)
    @testset "exp for Kronecker sum" begin
        EKS = exp(A ⊕ B)
        @test EKS isa AbstractKroneckerProduct
        @test EKS ≈ exp(kron(A, IB) + kron(IA, B))
    end

    @testset "sum over Kronecker sum" begin
        @test sum(KS) ≈ sum(kronsum)
        @test sum(KS3) ≈ sum(kronsum3)

        for dims in 1:2
            @test sum(KS, dims=dims) ≈ sum(kronsum, dims=dims)
            @test sum(KS3, dims=dims) ≈ sum(kronsum3, dims=dims)
        end

        @test_throws ArgumentError sum(KS, dims=-1)
        @test_throws ArgumentError sum(KS3, dims=3)
    end

end


text:
module RDatasets
    using Reexport, RData, CSV, CodecZlib
    using Printf: @sprintf
    @reexport using DataFrames

    export dataset

    global __packages = nothing
    global __datasets = nothing

    include("dataset.jl")
    include("datasets.jl")
    include("packages.jl")
end

max_stars_repo_name: SCXFPROTOCOL/GameRewards-System

max_issues_repo_name: SCXFPROTOCOL/GameRewards-System

max_forks_repo_name: SCXFPROTOCOL/GameRewards-System


text:
macro load(algo)
    esc(_load(string(algo)))
end
macro load(algos...)
    exprs = map(algos) do algo
        :(Nonconvex.@load $algo)
    end
    return Expr(:block, exprs...)
end

function _load(algo)
    if algo in ("MMA", "GCMMA", "MMA87", "MMA02")
        return install_and_load_module(:NonconvexMMA)
    elseif algo in ("Ipopt", "IpoptAlg")
        return install_and_load_module(:NonconvexIpopt)
    elseif algo in ("NLopt", "NLoptAlg")
        return install_and_load_module(:NonconvexNLopt)
    elseif algo in ("Juniper", "JuniperIpopt", "JuniperIpoptAlg")
        return install_and_load_module(:NonconvexJuniper)
    elseif algo in ("Pavito", "PavitoIpoptCbc", "PavitoIpoptCbcAlg")
        return install_and_load_module(:NonconvexPavito)
    elseif algo in ("Percival", "AugLag", "PercivalAlg")
        return install_and_load_module(:NonconvexPercival)
    elseif algo == "AugLag2"
        return install_and_load_module(:NonconvexAugLagLab)
    elseif algo in ("Bayesian", "BayesOpt", "BayesOptAlg")
        return install_and_load_module(:NonconvexBayesian)
    elseif algo in ("SDP", "Semidefinite", "SDPBarrier", "SDPBarrierAlg")
        return install_and_load_module(:NonconvexSemidefinite)
    elseif algo in ("Search", "MTS", "LS1", "MTSAlg", "LS1Alg")
        return install_and_load_module(:NonconvexSearch)
    elseif algo in ("Hyperopt", "Deflated", "Multistart", "HyperoptAlg", "DeflatedAlg")
        return install_and_load_module(:NonconvexMultistart)
    else
        throw("Unsupported algorithm. Please check the documentation of Nonconvex.jl.")
    end
end
function install_and_load_module(mod)
    quote
        using Pkg
        modname = $(Meta.quot(mod))
        try
            @info "Attempting to load the package $modname."
            using $mod
            @info "Loading succesful."
        catch
            @info "Couldn't find the package $modname. Attempting to install it."
            try
                Pkg.add(string(modname))
            catch
                @info "Package installation failed! Please report an issue."
                return
            end
            @info "$modname installed."
            @info "Attempting to load the package $modname."
            using $mod
            @info "Loading succesful."
        end
    end
end


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

mutable struct ImageDiskReference <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun

    function ImageDiskReference(;id=nothing, lun=nothing)
        o = new()
        validate_property(ImageDiskReference, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ImageDiskReference, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        o
    end
end # type ImageDiskReference

const _property_map_ImageDiskReference = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("lun")=>Symbol("lun"))
const _property_types_ImageDiskReference = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("lun")=>"Int32")
Base.propertynames(::Type{ ImageDiskReference }) = collect(keys(_property_map_ImageDiskReference))
Swagger.property_type(::Type{ ImageDiskReference }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ImageDiskReference[name]))}
Swagger.field_name(::Type{ ImageDiskReference }, property_name::Symbol) =  _property_map_ImageDiskReference[property_name]

function check_required(o::ImageDiskReference)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    true
end

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


text:
module Statistics

include("distributions.jl")

end


text: Base.show(io::IO, data::ModelSelection.ModelSelectionData) = ""


text:
function bfgs(nlp     :: AbstractNLPModel;
              x       :: Vector{T}=copy(nlp.meta.x0),
              ϵ       :: T = 1e-6,
              maxiter :: Int = 200,
              scaling :: Bool = true,
              Lp      :: Real = 2, # norm Lp
              B₀      :: Union{AbstractLinearOperator,
                               AbstractMatrix,
                               UniformScaling{T},
                               Nothing}              = nothing
              ) where T

    @info log_header([:iter, :f, :dual, :step, :slope], [Int, T, T, T, T],
                     hdr_override=Dict(:f=>"f(x)", :dual=>"‖∇f‖", :slope=>"∇fᵀd"))
    f = obj(nlp,x)
    ∇f = grad(nlp, x)

    xt = similar(x)
    ∇ft = similar(∇f)

    B = one(eltype(x))*I
    if B₀ != nothing
        B = B₀
    end
    B₀ = B
    #@show typeof(B)
    #@show B
    τ₀ = 0.0005
    τ₁ = 0.9999

    iter = 0
    @info log_row(Any[iter, f, norm(∇f)])

    while (norm(∇f, Lp) > ϵ) && (iter <= maxiter)
        d = - B*∇f

        #------------------------------------------
        # Hard coded line search
        hp0 = ∇f'*d
        t=1.0
        # Simple Wolfe forward tracking
        xt = x + t*d
        ∇ft = grad(nlp,xt)
        hp = ∇ft'*d
        ft = obj(nlp, xt)
        #  while  ~wolfe & armijo
        nbW = 0
        while (hp <= τ₁ * hp0) && (ft <= ( f + τ₀*t*hp0)) && (nbW < 10)
            t *= 5
            xt = x + t*d
            ∇ft = grad(nlp,xt)
            hp = ∇ft'*d
            ft = obj(nlp, xt)
            nbW += 1
            @debug "W", ft
        end
        tw = t

        # Simple Armijo backtracking
        nbk = 0
        while (ft > ( f + τ₀*t*hp0)) && (nbk < 20)
            t *= 0.4
            xt = x + t*d
            ft = obj(nlp, xt)
            nbk += 1
            @debug "A", ft
        end
        #------------------------------------------

        if t!=tw   ∇ft = grad(nlp, xt) end

        # Update BFGS approximation.
        B = push!(B, t * d, ∇ft - ∇f)
        #B = push!(B, t * d, ∇ft - ∇f, scaling)

        #move on
        x .= xt
        f = ft
        ∇f .= ∇ft

        iter += 1

        @info log_row(Any[iter, ft, norm(∇ft, Lp), t, hp0])
    end

    if iter > maxiter
        @warn "Maximum d'itérations"
        #@show B.data.M
    end

    return iter, f, norm(∇f, Lp), B, x
end

""" L-BFGS wrapper of the gereral BFGS implementation.
"""
function L_bfgs(nlp     :: AbstractNLPModel;
                x       :: Vector{T}=copy(nlp.meta.x0),
                ϵ       :: T = 1e-6,
                mem     :: Int = 5,
                maxiter :: Int = 200,
                scaling :: Bool = true,
                Lp      :: Real = 2, # norm Lp
                kwargs...
                ) where T

    @debug "U_Solver = L_bfgs"
    n = nlp.meta.nvar
    B₀ =  InverseLBFGSOperator(Float64, n, mem=mem, scaling=scaling)

    return bfgs(nlp; x=x, ϵ = ϵ, maxiter = maxiter, scaling=scaling, Lp = Lp, B₀=B₀, kwargs...)
end

""" M-BFGS wrapper of the gereral BFGS implementation.  M == Matrix
"""
function M_bfgs(nlp     :: AbstractNLPModel;
                x       :: Vector{T}=copy(nlp.meta.x0),
                ϵ       :: T = 1e-6,
                maxiter :: Int = 200,
                scaling :: Bool = true,
                Lp      :: Real = 2, # norm Lp
                kwargs...
                ) where T

    @debug "U_Solver = M_bfgs"
    n = nlp.meta.nvar
    B₀ =  InverseBFGSOperator(Float64, n, scaling=scaling)

    return bfgs(nlp; x=x, ϵ = ϵ, maxiter = maxiter, scaling=scaling, Lp = Lp, B₀=B₀, kwargs...)
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
module KnetOnnx

using ProtoBuf, MacroTools, Statistics, DataFlow, Knet

include("graph/onnx_pb.jl")
include("graph/new_types.jl")
include("graph/graph.jl")
include("graph/convert.jl")

include("converters.jl"); export ONNXtoGraph, PrintGraph;
include("KnetModel.jl"); export KnetModel;

end


text:
# Performance and runtime benchmark to the original MATLAB implementation

# julia> versioninfo()
# Julia Version 1.7.0-beta4
# Commit d0c90f37ba (2021-08-24 12:35 UTC)
# Platform Info:
#   OS: macOS (x86_64-apple-darwin19.6.0)
#   CPU: Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz
#   WORD_SIZE: 64
#   LIBM: libopenlibm
#   LLVM: libLLVM-12.0.1 (ORCJIT, skylake)
# Environment:
#   JULIA_NUM_THREADS = 8

using WNNMDenoise
using MAT, FileIO
using ImageQualityIndexes

using LinearAlgebra
using MKL
@info BLAS.get_config() # MKL is faster on svd

img = Float32.(load("house.png")) * 255;

σₙ = 40
noisy_img = matopen("house_AWGN_$(σₙ).mat") do io
    read(io, "N_Img")
end

assess_psnr(noisy_img, img, 255)
# MATLAB: 16.06
# Julia: 16.06

out = copy(noisy_img)
# f = WNNM(σₙ, K=2)
f = WNNM(σₙ)
# @time f(out, noisy_img; clean_img=img);
@time f(out, noisy_img);
# MATLAB: 97.365962 seconds
# Julia:  17.232886 seconds

assess_psnr(out, img, 255)
# MATLAB: 31.31
# Julia:  31.35


text:
module SecretSanta

function secret_santa(people::Array{String}, path::String = ENV["HOME"])
    n_people = length(people)

    if n_people < 2
        error("Need to provide at least two people to exchange gifts")
    end

    srand()
    perm = randperm(n_people)
    while some_same(perm) # Make sure no-one gifts themselves
        perm = randperm(n_people)
    end

    for i ∈ 1:n_people
        giver = people[i]
        receiver = people[perm[i]]

        open(joinpath(path, "$(giver).txt"), "w") do io
            write(io, """
                Dear $(giver),

                For Secret Santa, please provide a gift for $(receiver).

                Kind regards,
                Julia
                """)
        end
    end
end

function some_same(perm)
    for i ∈ 1:length(perm)
        if i == perm[i]
            return true
        end
    end
    return false
end

end # module


text:
struct SccNodeData
    depth::Int
    lowlink::Int
    onstack::Bool
end

struct StronglyConnectedComponent
    members::Vector{Int}
    is_bottom::Bool
end

export strongly_connected_components, StronglyConnectedComponent

function update_link(node_info::Vector{SccNodeData}, node::Int, new_node::Int)
    if node_info[node].lowlink > node_info[new_node].lowlink
        node_info[node] = SccNodeData(node_info[node].depth, node_info[new_node].lowlink, node_info[node].onstack)
    end
end

function strongly_connected_components(dg::Digraph)::Vector{StronglyConnectedComponent}
    components = Vector{StronglyConnectedComponent}()
    node_info = fill(SccNodeData(0, 0, false), node_count(dg))
    node_stack = Vector{Int}()
    depth = Int(1)
    function scc_recur(src::Int)
        node_info[src] = SccNodeData(depth, depth, true)
        depth += 1
        push!(node_stack, src)
        is_bottom = true
        for edge ∈ out_edges(dg, src)
            dst = dst_node(dg, edge)
            if  dst != src
                if node_info[dst].depth == 0
                    if !scc_recur(dst)
                        is_bottom = false
                    end
                    update_link(node_info, src, dst)
                elseif node_info[dst].onstack
                    update_link(node_info, src, dst)
                else
                    is_bottom = false
                end
            end
        end
        if node_info[src].lowlink == node_info[src].depth
            mem = Vector{Int}()
            local top_node::Int
            while true
                top_node = pop!(node_stack)
                push!(mem, top_node)
                node_info[top_node] = SccNodeData(node_info[top_node].depth, node_info[top_node].lowlink, false)
                if top_node == src
                    break
                end
            end
            push!(components, StronglyConnectedComponent(mem, is_bottom))
            return false
        else
            return is_bottom
        end
    end
    for n in nodes(dg)
        if node_info[n].depth == 0
            scc_recur(n)
        end
    end
    return components
end


text:
using CSV
using DataFrames
using Glob
using HDF5

import Base.@kwdef

## These should maybe go in a types file
"""
	PetaloDF

Returns Petalo data sets as DataFrames
"""
struct PetaloDF
	volume_names::Vector{String}
	process_names::Vector{String}
	sensor_xyz::DataFrame
	primaries::DataFrame
	vertices::DataFrame
	total_charge::DataFrame
	waveform::DataFrame
end

"""
    Define type for output so general dtype generator can be used easily.
"""
abstract type OutputDataset end

"""
	MlemLor

Struct representing a LOR
"""
struct MlemLor <: OutputDataset
    dt::Float32
    x1::Float32
    y1::Float32
    z1::Float32
    x2::Float32
    y2::Float32
    z2::Float32
	q1::Float32
	q2::Float32
	E1::Float32
	E2::Float32
end

"""
    EventParameters

Struct representing the parameters for each hemisphere of a selected event.
@kwdef macro allows for keyword argument definitions, improving readability.
"""
@kwdef struct EventParameters <: OutputDataset
    event_id ::Int64   = zero(Int64)
    phot1    ::Bool    = false
    phot2    ::Bool    = false
    nsipm1   ::Int64   = zero(Int64)
    nsipm2   ::Int64   = zero(Int64)
	q1       ::Float32 = zero(Float32)
    q2       ::Float32 = zero(Float32)
	E1       ::Float32 = zero(Float32)
	E2       ::Float32 = zero(Float32)
	r1       ::Float32 = zero(Float32)
    r2       ::Float32 = zero(Float32)
    r1x      ::Float32 = zero(Float32)
    r2x      ::Float32 = zero(Float32)
	phistd1  ::Float32 = zero(Float32)
    phistd2  ::Float32 = zero(Float32)
	widphi1  ::Float32 = zero(Float32)
    widphi2  ::Float32 = zero(Float32)
	zstd1    ::Float32 = zero(Float32)
    zstd2    ::Float32 = zero(Float32)
	widz1    ::Float32 = zero(Float32)
    widz2    ::Float32 = zero(Float32)
	corrzphi1::Float32 = zero(Float32)
    corrzphi2::Float32 = zero(Float32)
	xt1      ::Float32 = zero(Float32)
    xt2      ::Float32 = zero(Float32)
	yt1      ::Float32 = zero(Float32)
    yt2      ::Float32 = zero(Float32)
	zt1      ::Float32 = zero(Float32)
    zt2      ::Float32 = zero(Float32)
	ti1      ::Float32 = zero(Float32)
	ti2      ::Float32 = zero(Float32)
	t1       ::Float32 = zero(Float32)
    t2       ::Float32 = zero(Float32)
	x1       ::Float32 = zero(Float32)
    x2       ::Float32 = zero(Float32)
	y1       ::Float32 = zero(Float32)
    y2       ::Float32 = zero(Float32)
	z1       ::Float32 = zero(Float32)
    z2       ::Float32 = zero(Float32)
	xr1      ::Float32 = zero(Float32)
    xr2      ::Float32 = zero(Float32)
	yr1      ::Float32 = zero(Float32)
    yr2      ::Float32 = zero(Float32)
	zr1      ::Float32 = zero(Float32)
    zr2      ::Float32 = zero(Float32)
	tr1      ::Float32 = zero(Float32)
    tr2      ::Float32 = zero(Float32)
	xb1      ::Float32 = zero(Float32)
    xb2      ::Float32 = zero(Float32)
	yb1      ::Float32 = zero(Float32)
    yb2      ::Float32 = zero(Float32)
	zb1      ::Float32 = zero(Float32)
    zb2      ::Float32 = zero(Float32)
	ta1      ::Float32 = zero(Float32)
    ta2      ::Float32 = zero(Float32)
	xs       ::Float32 = zero(Float32)
    ux       ::Float32 = zero(Float32)
	ys       ::Float32 = zero(Float32)
    uy       ::Float32 = zero(Float32)
	zs       ::Float32 = zero(Float32)
    uz       ::Float32 = zero(Float32)
end

"""
    SimConfiguration

Type to store information about the simulation and analysis done on the data.
Will be extended when full electronics simulation ready!
"""
@kwdef struct SimConfiguration <: OutputDataset
    SiPMPDE::Float32 =    0.3
    ElecSTD::Float32 =    0.05
    SiPMThr::Float32 =    2.0
    SumQmin::Float32 =  100.0
    SumQmax::Float32 = 5000.0
    Ntof   ::Int64   =    7
    NEvent ::Int64   =    0
    Rmin   ::Float32 =  353.0
    Rmax   ::Float32 =  393.0
    CalLin ::Float32 =    0.0 #DOI Calibration linear term, to be implemented
    CalCon ::Float32 =    0.0 #DOI Calibration constant term
    CalBias::Float32 =    0.0 #DOI Calibration bias.
end

"""
	EventTypes

Type to save the counts of types of events in a dataset.
Contains:
total: Number of events generated.
empty: Number of events without any interaction.
single: Number of events with only one interaction.
prompt: Number of events with 2 reconstructable interactions.
single_prompt: Number of events with one interaction with q in selected range.
good_prompt: Number of events with two interactions with q in selected range.
"""
@kwdef mutable struct EventTypes <: OutputDataset
	total        ::Int64 = 0
	empty        ::Int64 = 0
	single       ::Int64 = 0
	prompt       ::Int64 = 0
	single_prompt::Int64 = 0
	good_prompt  ::Int64 = 0
end

## Readers/Writers
"""
    readdf(dir)

Reads all csv files found in dir and returns a single DataFrame.
"""
function readdf(dir::String)
    file_names = glob("*.csv", dir)
    dfs        = [CSV.read(fn, DataFrame) for fn in file_names]
    vcat(dfs...)
end

"""
    writemdf(dir::String, file::String, df::DataFrame)
    #NAME!!
writes a DataFrame to csv file.
"""
function writemdf(dir::String, file_name::String, df::DataFrame)
    path = joinpath(dir, file_name)
    CSV.write(path, df)
end

"""
	readh5_dset(path::String, folder::String, dset::String)

read an h5 dataset
"""
function readh5_dset(path::String, folder::String, dset::String)
	h5open(path, "r") do fid
		readh5_dset(fid, folder, dset)
	end
end

"""
	readh5_dset(h5in::HDF5.File, folder::String, dset::String)

read dataset from open hdf5 file.
"""
function readh5_dset(h5in::HDF5.File, folder::String, dset::String)
	g   = h5in[folder]
	dst = read(g, dset)
	return dst
end

"""
    readh5_todf(path::String, folder::String, dset::String)
Read contents of folder/dset in file path into a DataFrame.
"""
function readh5_todf(path::String, folder::String, dset::String)
    h5open(path, "r") do fid
		readh5_todf(fid, folder, dset)
	end
end

"""
    readh5_todf(h5in::HDF5.File, folder::String, dset::String)
Read from an open HDF5 file the data in folder/dset into a DataFrame.
"""
function readh5_todf(h5in::HDF5.File, folder::String, dset::String)
    df = readh5_dset(h5in, folder, dset) |> DataFrame
    return df
end

"""
	read_abc(path::String)

read abracadabra hdf5 and return relevant data sets in a PetaloDF struct.
"""
function read_abc(path::String)
	(primaries, sensor_xyz,
	 total_charge, vertices,
	 waveform, volume_names, process_names) = h5open(path, "r") do h5in
		a = readh5_todf(h5in, "MC", "primaries"    )
		b = readh5_todf(h5in, "MC", "sensor_xyz"   )
		c = readh5_todf(h5in, "MC", "total_charge" )
		d = readh5_todf(h5in, "MC", "vertices"     )
		e = readh5_todf(h5in, "MC", "waveform"     )
		f = readh5_dset(h5in, "MC", "volume_names" )
		g = readh5_dset(h5in, "MC", "process_names")
		return a, b, c, d, e, f, g
	end

	## Are these conversions really necessary?
	function int_conv(cols::Vector{Symbol})
		cols .=> ByRow(Int)
	end

	transform!(primaries , int_conv([:event_id ]), renamecols=false)
	transform!(sensor_xyz, int_conv([:sensor_id]), renamecols=false)
	transform!(total_charge,
		int_conv([:event_id, :sensor_id, :charge]), renamecols=false)
	transform!(waveform,
		int_conv([:event_id, :sensor_id]), renamecols=false)
	transform!(vertices,
		int_conv([:event_id, :track_id, :parent_id, :process_id, :volume_id]),
		renamecols=false)

	return PetaloDF(volume_names,
	                process_names,
					sensor_xyz,
					primaries,
					vertices,
					total_charge,
					waveform)
end

"""
    read_evtpar(infiles::Vector{String})

Get the selected event parameters from file along with normalisation
and rmin rmax info.
"""
function read_evtpar(infiles::Vector{String})
    rmin   = typemax(Float64)
    rmax   = typemin(Float64)
    norm   = zero(Int64)
    df_vec = DataFrame[]
    for fn in infiles
        s_conf = h5open(fn) do h5in
            sim_conf = readh5_dset(h5in, "configuration", "RunConfiguration")[1]
            push!(df_vec, readh5_todf(h5in, "selected_events", "EventParameters"))
            return sim_conf
        end
        norm += s_conf.NEvent
        rmin  = s_conf.Rmin < rmin ? s_conf.Rmin : rmin
        rmax  = s_conf.Rmax > rmax ? s_conf.Rmax : rmax
    end
    return norm, rmin, rmax, vcat(df_vec...)
end

"""
    generate_hdf5_datatype

returns a datatype object for saving any of the OutputDataset group type to hdf5.
"""
function generate_hdf5_datatype(data_type::Type{<:OutputDataset})
    epdtype = HDF5.h5t_create(HDF5.H5T_COMPOUND, sizeof(data_type))
    epnames = fieldnames(data_type)
    eptypes = fieldtypes(data_type)
    for (i, (nm, typ)) in enumerate(zip(epnames, eptypes))
        HDF5.h5t_insert(epdtype, nm, fieldoffset(data_type, i), datatype(typ))
    end
    return HDF5.Datatype(epdtype)
end

"""
	write_lors_hdf5(filename, mlor)
	Write lors in a hdf5 format required by petalorust (mlem algo)
"""
function write_lors_hdf5(filename, mlor, directory::String="reco_info")
	h5open(filename, "w") do h5f
		dtype  = generate_hdf5_datatype(MlemLor)
		dspace = dataspace(mlor)
		grp    = create_group(h5f, directory)
		dset   = create_dataset(grp, "lors", dtype, dspace)
		write_dataset(dset, dtype, mlor)
	end
end

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
# Slice sampler (Neal 2003; MacKay 2003, sec. 29.7)
#
# A slice sampler is an adaptive step-size MCMC algorithm for continuous random
# variables, which only requires an unnormalized density function as input.  It
# is a convenient alternative to Metropolis because it often gives good results
# with very little tuning, and works in the case that different parts of the
# distribution require different step sizes.  However, if the distribution's
# widths don't vary too much, and there's a good initialization and proposal,
# Metropolis may be more efficient.
#
# This slice sampler is univariate; for a multivariate problem it just updates
# variables one at a time.  Thus it suffers when variables are correlated.
#
# REQUIRED
#   logdist: Log-density function of target distribution
#   initial: Initial state
#   niter:   Number of iterations (samples to return)
#
# OPTIONAL
#   widths: Step sizes for initially expanding the slice (D-dim vector).
#           There is little harm in making these very large.
#   burnin: Set to >0 to run for 'burnin' iterations without recording values
#   step_out: Protects against the case if you pass in widths that are too small.
#             If you are sure your widths are large enough, can set this to false.
#   verbose: Set to true for information at every iteration
#
# RETURNS a sampling history.  Pick the last one for an independent sample.
#
# EXAMPLES: take samples from N(3,1), starting from a poor initializer. Mixes to a good place!
# Get many samples:
#   slice_sample(x-> -0.5(x-3)^2, 42.0, 30)
# Get one sample:
#   slice_sample(x-> -0.5(x-3)^2, 42.0, 30)[end]
#
# There are two versions of this procedure: either
#     univariate:   initial \in R,   logdist: R -> R,   retval \in R^niter
#     multivariate: initial \in R^D, logdist: R^D -> R, retval \in R^(niter x D)
#
# This is a port of Iain Murray's
# http://homepages.inf.ed.ac.uk/imurray2/teaching/09mlss/slice_sample.m

# (This code is the multivariate version; univariate interfaces are further below.)
function slice_sample(logdist::Function, initial::Array{Float64,1}, niter::Integer;
        widths::Array{Float64,1} = ones(length(initial)),
        step_out = true,
        burnin = 0,
        verbose = false)
    D = length(initial)
    @assert length(widths)==D

    state::Vector{Float64} = copy(initial)
    log_Px::Float64 = logdist(state)

    history = zeros(niter,D)

    for iter=1:(niter+burnin)
        if verbose
            @printf("Slice iter %d state [%s] log_Px %f\n",iter, join(map(x->@sprintf("%.5f",x),state)," "), log_Px)
        end

        # Sweep through axes
        for dd=1:D
            log_uprime = log(rand()) + log_Px
            x_l  = copy(state)
            x_r  = copy(state)
            xprime = copy(state)
            # Create a horizontal interval (x_l, x_r) enclosing xx
            r = rand()
            x_l[dd] = state[dd] - r*widths[dd]
            x_r[dd] = state[dd] + (1-r)*widths[dd]
            if step_out
                while logdist(x_l) > log_uprime
                    x_l[dd] -= widths[dd]
                end
                while logdist(x_r) > log_uprime
                    x_r[dd] += widths[dd]
                end
            end
            # Inner loop:
            # Propose xprimes and shrink interval until good one is found.
            while true
                xprime[dd] = rand() * (x_r[dd] - x_l[dd]) + x_l[dd];
                log_Px = logdist(xprime)
                if log_Px > log_uprime
                    break
                else
                    if (xprime[dd] > state[dd])
                        x_r[dd] = xprime[dd];
                    elseif (xprime[dd] < state[dd])
                        x_l[dd] = xprime[dd];
                    else
                        @assert false "BUG, shrunk to current position and still not acceptable";
                    end
                end
            end
            state[dd] = xprime[dd]
        end
        if iter > burnin
            history[iter-burnin, :] = state
        end
    end
    history
end

## Univariate formulations

function slice_sample(logdist::Function, initial::Float64, niter::Int; kwargs...)
    history = slice_sample(x-> logdist(x[1]), [initial], niter; kwargs...)
    reshape(history, (size(history,1),))
end


text:
# Raw memory management

export Mem, available_memory, total_memory

module Mem

using ..VectorEngine
using ..VectorEngine.VEDA: vedaMemAlloc, vedaMemPtr, vedaMemFree, vedaMemGetInfo,
    vedaMemAllocHost, vedaMemFreeHost

using Printf

#
# buffers
#

# a chunk of memory allocated using the VEDA APIs. this memory can reside on the host, on
# the VE, or can represent specially-formatted memory (like texture arrays). depending on
# all that, the buffer may be `convert`ed to a Ptr, VEPtr, or VEArrayPtr.

abstract type AbstractBuffer end

Base.convert(T::Type{<:Union{Ptr,VEPtr,VEArrayPtr}}, buf::AbstractBuffer) =
    throw(ArgumentError("Illegal conversion of a $(typeof(buf)) to a $T"))

# ccall integration
#
# taking the pointer of a buffer means returning the underlying pointer,
# and not the pointer of the buffer object itself.
Base.unsafe_convert(T::Type{<:Union{Ptr,VEPtr,VEArrayPtr}}, buf::AbstractBuffer) =
    convert(T, buf)

## host side device buffer

"""
    Mem.DeviceBuffer

Host residing structure representing a buffer of device memory.
"""
mutable struct DeviceBuffer <: AbstractBuffer
    #vptr::VEDAdeviceptr
    vptr::VEPtr{Int8}
    ptr::VEPtr{Int8}
    bytesize::Int
    function DeviceBuffer(bsize::Int)
        bsize == 0 && return new(VE_NULL, VE_NULL, 0)

        #vp = Ref{VEDAdeviceptr}()
        vp = Ref{VEPtr{Int8}}(0)
        p = Ref{VEPtr{Int8}}(0)
        vedaMemAlloc(pointer_from_objref(vp), bsize)
        vedaMemPtr(pointer_from_objref(p), vp[])
        #VectorEngine.vesync()      # is this needed?
        obj = new(vp[], p[], bsize)
        finalizer(unsafe_free!, obj)
        return obj
    end
end

function unsafe_free!(buf::DeviceBuffer)
    if pointer(buf) != VE_NULL
        vedaMemFree(pointer(buf))
        buf.vptr = VE_NULL
        buf.ptr = VE_NULL
    end
end

"""
    Mem.alloc(DeviceBuffer, bytesize::Integer)

Allocate `bytesize` bytes of memory on the VE.
"""
alloc(::Type{DeviceBuffer}, bytesize::Integer) = DeviceBuffer(bytesize)

function free(buf::DeviceBuffer) end

Base.pointer(buf::DeviceBuffer) = buf.vptr
Base.sizeof(buf::DeviceBuffer) = buf.bytesize

Base.show(io::IO, buf::DeviceBuffer) =
    @printf(io, "DeviceBuffer(%s at %p (%p))", Base.format_bytes(sizeof(buf)), Int(buf.vptr), Int(buf.ptr))

#Base.convert(::Type{VEDAdeviceptr}, buf::DeviceBuffer) = buf.vptr
Base.convert(::Type{VEPtr{T}}, buf::DeviceBuffer) where {T} = reinterpret(VEPtr{T}, buf.vptr)

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

## host buffer

"""
    Mem.HostBuffer
    Mem.Host

A buffer of pinned memory on the CPU, unaccessible to the VE.
"""
mutable struct HostBuffer <: AbstractBuffer
    ptr::Ptr{Cvoid}
    bytesize::Int
end

Base.pointer(buf::HostBuffer) = buf.ptr
Base.sizeof(buf::HostBuffer) = buf.bytesize

Base.show(io::IO, buf::HostBuffer) =
    @printf(io, "HostBuffer(%s at %p)", Base.format_bytes(sizeof(buf)), Int(pointer(buf)))

Base.convert(::Type{Ptr{T}}, buf::HostBuffer) where {T} =
    convert(Ptr{T}, pointer(buf))

function Base.convert(::Type{VEPtr{T}}, buf::HostBuffer) where {T}
    throw(ArgumentError("cannot take the VE address of a CPU buffer"))
end

"""
    Mem.alloc(HostBuffer, bytesize::Integer)

Allocate `bytesize` bytes of page-locked memory on the host. This memory is accessible from
the CPU, and makes it possible to perform faster memory copies to the VE.
"""
function alloc(::Type{HostBuffer}, bytesize::Integer)
    bytesize == 0 && return HostBuffer(C_NULL, 0)

    ptr_ref = Ref{Ptr{Cvoid}}()
    vedaMemAllocHost(ptr_ref, bytesize)

    return HostBuffer(ptr_ref[], bytesize)
end

function free(buf::HostBuffer)
    if pointer(buf) != VE_NULL
        vedaMemFreeHost(buf.ptr)
    end
end

## array buffer

mutable struct ArrayBuffer{T,N} <: AbstractBuffer
    ptr::VEArrayPtr{T}
    dims::Dims{N}
end

Base.pointer(buf::ArrayBuffer) = buf.ptr
Base.sizeof(buf::ArrayBuffer) = error("Opaque array buffers do not have a definite size")
Base.size(buf::ArrayBuffer) = buf.dims
Base.length(buf::ArrayBuffer) = prod(buf.dims)
Base.ndims(buf::ArrayBuffer{<:Any,N}) where {N} = N

Base.show(io::IO, buf::ArrayBuffer{T,1}) where {T} =
    @printf(io, "%g-element ArrayBuffer{%s,%g}(%p)", length(buf), string(T), 1, Int(pointer(buf)))
Base.show(io::IO, buf::ArrayBuffer{T}) where {T} =
    @printf(io, "%s ArrayBuffer{%s,%g}(%p)", Base.inds2string(size(buf)), string(T), ndims(buf), Int(pointer(buf)))

# array buffers are typed, so refuse arbitrary conversions
Base.convert(::Type{VEArrayPtr{T}}, buf::ArrayBuffer{T}) where {T} =
    convert(VEArrayPtr{T}, pointer(buf))
# ... except for VEArrayPtr{Nothing}, which is used to call untyped API functions
Base.convert(::Type{VEArrayPtr{Nothing}}, buf::ArrayBuffer)  =
    convert(VEArrayPtr{Nothing}, pointer(buf))

function alloc(::Type{<:ArrayBuffer{T}}, dims::Dims{N}) where {T,N}
    format = convert(CUarray_format, eltype(T))

    if N == 2
        width, height = dims
        depth = 0
        @assert 1 <= width "VEDA 2D array (texture) width must be >= 1"
        # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_WIDTH
        @assert 1 <= height "VEDA 2D array (texture) height must be >= 1"
        # @assert height <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_HEIGHT
    elseif N == 3
        width, height, depth = dims
        @assert 1 <= width "VEDA 3D array (texture) width must be >= 1"
        # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_WIDTH
        @assert 1 <= height "VEDA 3D array (texture) height must be >= 1"
        # @assert height <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_HEIGHT
        @assert 1 <= depth "VEDA 3D array (texture) depth must be >= 1"
        # @assert depth <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_DEPTH
    elseif N == 1
        width = dims[1]
        height = depth = 0
        @assert 1 <= width "VEDA 1D array (texture) width must be >= 1"
        # @assert witdh <= CU_DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE1D_WIDTH
    else
        "VEDA arrays (texture memory) can only have 1, 2 or 3 dimensions"
    end

    allocateArray_ref = Ref(CUDA.CUDA_ARRAY3D_DESCRIPTOR(
        width, # Width::Csize_t
        height, # Height::Csize_t
        depth, # Depth::Csize_t
        format, # Format::CUarray_format
        UInt32(CUDA.nchans(T)), # NumChannels::UInt32
        0))

    handle_ref = Ref{CUarray}()
    vedaArray3DCreate(handle_ref, allocateArray_ref)
    ptr = reinterpret(VEArrayPtr{T}, handle_ref[])

    return ArrayBuffer{T,N}(ptr, dims)
end

function free(buf::ArrayBuffer)
    vedaArrayDestroy(buf.ptr)
end

## convenience aliases

const Device  = DeviceBuffer
const Host    = HostBuffer
const Array   = ArrayBuffer

#
# pointers
#

## initialization

"""
    Mem.set!(buf::VEPtr, value::Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,
                                      Int64,UInt64,Float32,Float64}, len::Integer;
             async::Bool=false, stream::VEStream)

Initialize device memory by copying `val` for `len` times. Executed asynchronously if
`async` is true, in which case a valid `stream` is required.
"""
set!

for T in [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64]
    bits = 8*sizeof(T)
    fn_sync = Symbol("vedaMemsetD$(bits)")
    fn_async = Symbol("vedaMemsetD$(bits)Async")
    U = Symbol("UInt$(bits)")
    @eval function set!(ptr::VEPtr{$T}, value::$T, len::Integer;
                        async::Bool=false, stream::Union{Nothing,VEStream}=nothing)
        val = typeof(value) == $U ? value : reinterpret($U, value)
        if async
          stream===nothing &&
              throw(ArgumentError("Asynchronous memory operations require a stream."))
            $(getproperty(VectorEngine.VEDA, fn_async))(ptr, val, len, stream)
        else
          stream===nothing ||
              throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
            $(getproperty(VectorEngine.VEDA, fn_sync))(ptr, val, len)
        end
    end
end

## copy operations

for (f, fa, srcPtrTy, dstPtrTy) in (("vedaMemcpyDtoH", "vedaMemcpyDtoHAsync", VEPtr, Ptr),
                                    ("vedaMemcpyHtoD", "vedaMemcpyHtoDAsync", Ptr,   VEPtr),
                                    ("vedaMemcpyDtoD", "vedaMemcpyDtoDAsync", VEPtr, VEPtr),
                                   )
    @eval function Base.unsafe_copyto!(dst::$dstPtrTy{T}, src::$srcPtrTy{T}, N::Integer;
                                       stream::Union{Nothing,VEStream}=nothing,
                                       async::Bool=false) where T
        if async
            stream===nothing &&
                throw(ArgumentError("Asynchronous memory operations require a stream."))
            $(getproperty(VectorEngine.VEDA, Symbol(fa)))(dst, src, N*sizeof(T), stream)
        else
            stream===nothing ||
                throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
            $(getproperty(VectorEngine.VEDA, Symbol(f)))(dst, src, N*sizeof(T))
        end
        return dst
    end
end

function Base.unsafe_copyto!(dst::VEArrayPtr{T}, src::Ptr{T}, N::Integer;
                             stream::Union{Nothing,VEStream}=nothing,
                             async::Bool=false) where T
    if async
        stream===nothing &&
            throw(ArgumentError("Asynchronous memory operations require a stream."))
        vedaMemcpyHtoDAsync(dst, src, N*sizeof(T), stream)
    else
        stream===nothing ||
            throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
        vedaMemcpyHtoD(dst+doffs-1, src, N*sizeof(T))
    end
end

function Base.unsafe_copyto!(dst::Ptr{T}, src::VEArrayPtr{T}, soffs::Integer, N::Integer;
                             stream::Union{Nothing,VEStream}=nothing,
                             async::Bool=false) where T
    if async
        stream===nothing &&
            throw(ArgumentError("Asynchronous memory operations require a stream."))
        vedaMemcpyAtoHAsync(dst, src, soffs, N*sizeof(T), stream)
    else
        stream===nothing ||
            throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
        vedaMemcpyAtoH(dst, src, soffs, N*sizeof(T))
    end
end

Base.unsafe_copyto!(dst::VEArrayPtr{T}, doffs::Integer, src::VEPtr{T}, N::Integer) where {T} =
    vedaMemcpyDtoA(dst, doffs, src, N*sizeof(T))

Base.unsafe_copyto!(dst::VEPtr{T}, src::VEArrayPtr{T}, soffs::Integer, N::Integer) where {T} =
    vedaMemcpyAtoD(dst, src, soffs, N*sizeof(T))

Base.unsafe_copyto!(dst::VEArrayPtr, src, N::Integer; kwargs...) =
    Base.unsafe_copyto!(dst, 0, src, N; kwargs...)

Base.unsafe_copyto!(dst, src::VEArrayPtr, N::Integer; kwargs...) =
    Base.unsafe_copyto!(dst, src, 0, N; kwargs...)

## memory info

function info()
    free_ref = Ref{Csize_t}()
    total_ref = Ref{Csize_t}()
    vedaMemGetInfo(free_ref, total_ref)
    return convert(Int, free_ref[]), convert(Int, total_ref[])
end

end # module Mem

"""
    available_memory()

Returns the available_memory amount of memory (in bytes), available for allocation by the CUDA context.
"""
available_memory() = Mem.info()[1]

"""
    total_memory()

Returns the total amount of memory (in bytes), available for allocation by the CUDA context.
"""
total_memory() = Mem.info()[2]

# memory operations

function unsafe_fill!(ptr::Union{Ptr{T},VEPtr{T}},
                      pattern::Union{Ptr{T},VEPtr{T}}, N::Integer) where T
    bytes = N*sizeof(T)
    bytes==0 && return
    Mem.set!(ptr, pattern, N)
end


text:
import Printf: ini_dec, fix_dec, ini_hex, ini_HEX

function show(io::IO, x::DoubleFloat{T}) where {T<:IEEEFloat}
    compact = get(io, :compact, true)
    if compact
        str = string(x.hi)
    else
        str = string(x)
    end
    print(io, str)
end

function show(io::IO, x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    re, im = reim(x)
    imstr = isfinite(im) ? "im" : "*im"
    compact = get(io, :compact, true)
    if compact
        str = string(HI(re), (signbit(x.im.hi) ? " - " : " + "), abs(HI(im)), imstr)
    else
        str = string(re, (signbit(HI(im)) ? " - " : " + "), abs(im), imstr)
    end
    print(io, str)
end

show(x::DoubleFloat{T}) where {T<:IEEEFloat} = show(Base.stdout, x)

show(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat} = show(Base.stdout, x)

function showtyped(io::IO, x::DoubleFloat{T}) where {T<:IEEEFloat}
    str = stringtyped(x)
    print(io, str)
end

function showtyped(io::IO, x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat}
    str = stringtyped(x)
    print(io, str)
end

showtyped(x::DoubleFloat{T}) where {T<:IEEEFloat} = showtyped(Base.stdout, x)

showtyped(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat} = showtyped(Base.stdout, x)

showall(x::DoubleFloat{T}) where {T<:IEEEFloat} = print(Base.stdout, string(x))

showall(x::Complex{DoubleFloat{T}}) where {T<:IEEEFloat} = print(Base.stdout, string(x))

if VERSION < v"1.1"
    fix_dec(x::Double64, n::Int) = fix_dec(Float64(x), n)
    fix_dec(x::Double32, n::Int) = fix_dec(Float64(x), n)
    fix_dec(x::Double16, n::Int) = fix_dec(Float32(x), n)
    ini_dec(x::Double64, n::Int) = ini_dec(Float64(x), n)
    ini_dec(x::Double32, n::Int) = ini_dec(Float64(x), n)
    ini_dec(x::Double16, n::Int) = ini_dec(Float32(x), n)
else
    fix_dec(x::Double64, n::Int, digits) = fix_dec(Float64(x), n, digits)
    fix_dec(x::Double32, n::Int, digits) = fix_dec(Float64(x), n, digits)
    fix_dec(x::Double16, n::Int, digits) = fix_dec(Float32(x), n, digits)
    ini_dec(x::Double64, n::Int, digits) = ini_dec(Float64(x), n, digits)
    ini_dec(x::Double32, n::Int, digits) = ini_dec(Float64(x), n, digits)
    ini_dec(x::Double16, n::Int, digits) = ini_dec(Float32(x), n, digits)
end


text:
"""
    augment(term) -> Vector{String}

Given a term, returns a list of terms which should be treated as synonyms.
Currently only supports agumenting (spaces or hyphens) with (spaces, no spaces).

## Example

```jldoctest
julia> KeywordSearch.augment("arctic wolf")
2-element Vector{String}:
 "arctic wolf"
 "arcticwolf"

```
"""
function augment(term)
    terms = String[]
    words = split(term, char -> isspace(char) || char == '-')
    n_words = length(words)
    joiners = (" ", "")
    # This could be optimized
    for joins in Iterators.product(Iterators.repeated(joiners, n_words - 1)...)
        word = ""
        for (w, j) in zip(words, joins)
            word = string(word, w, j)
        end
        word = string(word, words[end])
        push!(terms, word)
    end
    return terms
end

reconstruct(Q::FuzzyQuery) = x -> FuzzyQuery(x, Q.dist, Q.threshold)
reconstruct(Q::Query) = Query

function augment(Q::AbstractQuery)
    make_q = reconstruct(Q)
    return Or(Tuple(make_q.(augment(Q.text))))
end

# `raw` needed due to the `\n` in the test
@doc raw"""
    word_boundary(Q::AbstractQuery) -> AbstractQuery

Ensures that a word or phrase is not hyphenated or conjoined with the surrounding text.

## Example

```jldoctest
julia> using Test

julia> query = Query("word")
Query("word")

julia> @test match(query, Document("This matchesword ")) !== nothing
Test Passed

julia> @test match(word_boundary(query), Document("This matches word.")) !== nothing
Test Passed

julia> @test match(word_boundary(query), Document("This matches word ")) !== nothing
Test Passed

julia> @test match(word_boundary(query), Document("This matches word\nNext line")) !== nothing
Test Passed

julia> @test match(word_boundary(query), Document("This doesn't matchword ")) === nothing
Test Passed

```
"""
function word_boundary(Q::AbstractQuery)
    # `process_document` has removed punctuation, so we just need to check for spaces.
    make_q = reconstruct(Q)
    stripped_text = strip(Q.text)
    return make_q(string(" ", stripped_text, " "))
end

Base.length(Q::AbstractQuery) = mapreduce(x -> 1, +, AbstractTrees.Leaves(Q); init=0)

function check_keys(::Type{T}) where {T<:NamedTuple}
    if :match ∈ _nt_names(T)
        throw(ArgumentError("Must not include `match` as a metadata key."))
    end
    return nothing
end

function disjoint_keys_check(::NamedQuery{T1}, ::Corpus{T2,TR}) where {T1,T2,TR}
    disjoint_keys_check(T1, T2)
    disjoint_keys_check(T1, TR)
    return nothing
end

function disjoint_keys_check(::NamedQuery{T1}, ::Document{T2}) where {T1,T2}
    disjoint_keys_check(T1, T2)
    return nothing
end

# a similar version of this exists as an internal function in Base:
# https://github.com/JuliaLang/julia/blob/e68dda9785b4523cae49f8a60f99aa9360226eb4/base/namedtuple.jl#L181
_nt_names(::Type{NamedTuple{names,T}}) where {names,T} = names

function disjoint_keys_check(::Type{T1}, ::Type{T2}) where {T1,T2}
    if !isdisjoint(_nt_names(T1), _nt_names(T2))
        error("Metadata keys will clash when merging metadata to construct the `NamedMatch`.")
    end
    return nothing
end

max_stars_repo_name: beacon-biosignals/KeywordSearch

max_issues_repo_name: beacon-biosignals/KeywordSearch

max_forks_repo_name: beacon-biosignals/KeywordSearch


text:
using RegistryCompatTools

using Test

const rexversion = "\\d+\\.\\d+\\.\\d+"

@test held_back_packages() isa Dict
@test held_back_by("Requires") isa Vector{String}
let d = held_back_packages()
    io = IOBuffer()
    show(IOContext(io, :color=>true), MIME("text/plain"), first(values(d)))
    str = String(take!(io))
    pkgs = split(chomp(str), '\n')[2:end]
    @test !isempty(pkgs)
    @test occursin(Regex(".*@.*$rexversion.* .*$rexversion"), first(pkgs))   # check for version string (allow ANSI skip)
    @test occursin("\e[", first(pkgs))   # test ANSI terminal codes for coloration
end

# The following test assumes that at least one of the dependents of ColorTypes is compatible with the latest version
@test length(held_back_by("ColorTypes", typemax(VersionNumber))) > length(held_back_by("ColorTypes"))


text:
### A Pluto.jl notebook ###
# v0.12.17

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 935e57b0-438f-11eb-0ed4-b3ad455bbd63
using STMOZOO.BeesAlgorithm

# ╔═╡ 0ca837e0-42ef-11eb-17fa-9335cb9a3997
using InteractiveUtils, Plots, PlutoUI

# ╔═╡ 91d40c30-43a5-11eb-1805-ef48b11c3a40
using Optim

# ╔═╡ eee523f0-436e-11eb-2e92-d59ed5c533e6
md"
# The Artificial Bee Algorithm - Tutorial

## Introduction
The artificial bee algorithm (ABC) is a swarm intelligence metaheuristics algorithm for solving optimization problems. In this tutorial, the module **BeesAlgorithm**, which implements the ABC algorithm for optimization of **continuous functions**, will be demonstrated. Several test functions such as the **Sphere, Ackley, Rosenbrock, Branin** and **Rastrigine** function will be minimized to illustrate the efficacy of the ABC algorithm. Visualizations are provided to show the convergence of the ABC algorithm.

## Concept
##### General
The ABC algorithm is inspired by the **foraging behaviour of honey bees**. Honey bees collect nectar from flower patches as a food source for the colony. Bees are sent to explore different flower patches and they communicate the quality of these food sources through waggle dances. Good sites are continually exploited, while bees are sent out in search of additional promising sites.

##### Metaphor

As a methaphor for the foraging behaviour of bees, the ABC algorithm relies on 3 main components:

- **Food sources**, which can be considered as potential solutions of the optimization problem.

- **Employed foragers**. They exploit a food source, return to the colony to share their information with a certain probability, perform a waggle dance and recruit other bees, and then continue to forage at the food source.

- **Unemployed foragers**. This category consists of 2 types of bees. On the one hand, the **onlooker bees** watch the waggle dances to become a recruit and start searching for a food source. On the other hand, the **scout bees** start searching for interesting flower patches around the nest spontaneously.

The fitness of a solution or food source is inversely related with the value of the objective function in this solution. Thus, a higher fitness corresponds to a lower objective value. In the optimization process, we want to **maximize fitness** and **minimize the objective function** to find the minimizer of a continous function.

The following phases in the ABC algorithm can be distinguished:

**1) Employed bee phase**\
Employed bees try to identify better food source than the one they were associated previously. A new solution is generated using a partner solution. Thereafter, greedy selection is performed, meaning that a new solution only  will be accepted if it is better than the current solution. Every bee in the swarm will explore one food source. All solutions get an opportunity to generate a new solution in the employed bee phase.

**2) Onlooker bee phase**\
In the onlooker bee phase, a food source is selected for further exploitation with a probability related to the nectar amount, i.e. a solution with higher fitness will have a higher probability to be chosen. Fitter solutions may undergo multiple onlooker bee explorations. As in the employed bee phase, new solutions are generated using a partner solution and greedy selection is performed. In contrast to the employed bee phase, not every food source will be explored, since every onlooker bee will explore a certain food source with a certain probability (depending on nectar amount).

During the 2 phases above, a trial counter is registered for every food source. Each time a food source fails to generate a solution with higher fitness, the trial counter is elevated by 1 unit.

The solution with highest fitness so far is kept apart in memory during the entire process and updated as better food sources are discovered.

**3) Scout bee phase**\
If the value of the trial counter for a certain solution is greater than fixed limit, then a solution can enter the scout phase. The latter food source is then considered as exhausted and will therefore be abandoned by the bees. After discarding the exhausted solution, a new random solution is generated and the trial counter of this solution is reset to zero.

"

# ╔═╡ 3c1c9be0-4373-11eb-1b1e-3d3f2e6afb19
md"
## Visualization of the ABC algorithm
#### Evolution of bee swarms during optimization

Type *add https://github.com/kirstvh/STMOZOO.git* in the package manager of the Julia terminal to download the Bees Algorithm module.
"

# ╔═╡ 7a2fdff0-4373-11eb-337a-8fb102ef7b78
md"  **Loading the Bees Algorithm module** "

# ╔═╡ 16e1ad62-439c-11eb-1b3e-577f23587233
md"
The explanation of the functions in the STMOZOO.BeesAlgorithm module can be retrieved from the documentation page or by using the help function."

# ╔═╡ 02f4f820-439c-11eb-0c6e-7509dfc20674
md"  **Loading the required packages for this tutorial** "

# ╔═╡ a9462a5e-4373-11eb-3b48-39a2596229a9
md" **Choose the test function** you want to minimize with the select button below."

# ╔═╡ 27f302ee-42ea-11eb-2d9e-49dffc0d983d
@bind functie Select(["ackley", "sphere","rosenbrock","branin","rastrigine"])

# ╔═╡ 3235a8d2-42ea-11eb-1fe1-6d91eca83dad
begin
	if functie == "ackley"
		f_optimize = ackley;
	end
	if functie == "sphere"
		f_optimize = sphere;
	end
	if functie == "rosenbrock"
		f_optimize = rosenbrock;
	end
	if functie == "branin"
		f_optimize = branin;
	end
	if functie == "rastrigine"
		f_optimize = rastrigine;
	end
end

# ╔═╡ 085c34e0-4374-11eb-1ba7-7fb1af1d38a4
md" Choose the **swarm size S** (even number). \
 "

# ╔═╡ 0388b3ce-4374-11eb-03f4-9b4c74bd5ff0
S = 30

# ╔═╡ 6bd962e0-439c-11eb-0ce9-c348bd45b225
md"
The **number of decision variables D** (determining the dimension of the optimization problem) is fixed to 2 in this tutorial for  visualization purposes. The **limit** is fixed to a convenient value of D*(S/2). "

# ╔═╡ 350644d0-4375-11eb-2b68-8fb8f4bd7c2a
D = 2

# ╔═╡ c7c64720-437f-11eb-15a2-477ab2fe0792
limit = D * (S/2)

# ╔═╡ caa5b6d0-4373-11eb-0cdd-1961e6698727
md" Appropriate parameters for the chosen test function are defined. Next, the **ABC algorithm executes the searching procedure** for finding the global minimum of the test function."

# ╔═╡ f347e610-42a3-11eb-2116-ef50f1246cf3
begin

	if functie == "sphere"
		T = 35
		bounds_lower = [-100,-100];
		bounds_upper = [100,100];
	end

	if functie == "ackley"
		T = 50
		bounds_lower = [-30,-30];
		bounds_upper = [30,30];
	end

	if functie == "rosenbrock"
		T = 1000
		bounds_lower = [-100, -100];
		bounds_upper = [100,100];
	end

	if functie == "branin"
		T = 50
		bounds_lower = [-5,0];
		bounds_upper = [10,15];
	end

	if functie == "rastrigine"
		T = 50
		bounds_lower = [-5,-5];
		bounds_upper = [5,5];
	end
end

# ╔═╡ df912bb0-43a5-11eb-3bd2-ebc17c073757
begin
	using STMOZOO.BeesAlgorithm: initialize_population
	starting_solution = float(initialize_population(D, bounds_lower, bounds_upper, S/2)[1])
end

# ╔═╡ 54c02380-42a4-11eb-0240-7b2d895cb337
begin
	optimal_solution,populations, best_fitness_tracker = ArtificialBeeColonization(D, bounds_lower, bounds_upper, S, T, limit, f_optimize)
	optimal_solution
end

# ╔═╡ 9bb28f50-4374-11eb-2b10-e5effcbc8438
md" Below a **contour plot** and **surface plot** of the test function can be seen. The animation shows how the location of bees changes over time during the optimization procedure. "

# ╔═╡ b81d7f30-42a5-11eb-27ce-f1cc849ffdc5

# ╔═╡ ae391ef0-47d4-11eb-15ae-a7d5f188644b
begin

			if functie == "sphere"
			x2=range(bounds_lower[1],bounds_upper[1], step=1)
			y2=range(bounds_lower[2],bounds_upper[2], step=1)
			f(x2,y2) = (x2^2+y2^2)
		end

		if functie == "ackley"
			x2=range(bounds_lower[1],bounds_upper[1], step=0.75)
			y2=range(bounds_lower[2],bounds_upper[2], step=0.75)
		    d = 2
			c=2*3.14
			a=20
			b=0.2
		    f(x2,y2) = -a * exp(-b*sqrt((x2^2+y2^2)/d))-exp((cos(c*x2)+cos(c*y2))/d) + a + exp(1)

		end

		if functie == "rosenbrock"
			x2=range(bounds_lower[1],bounds_upper[1], step=0.5)
			y2=range(bounds_lower[2],bounds_upper[2], step=0.5)
			a=1
			b=5
		    f(x2,y2) = (a-x2)^2 + b*(y2-x2^2)^2
		end

		if functie == "branin"
			x2=range(bounds_lower[1],bounds_upper[1], step=0.5)
			y2=range(bounds_lower[2],bounds_upper[2], step=0.5)
			a=1
			b=5.1/(4pi^2)
			c=5/pi
			r=6
			s=10
			t=1/8pi
		    f(x2,y2) = a * (y2 - b * x2^2 + c * x2 - r)^2 + s * (1 - t) * cos(x2) + s
		end

		if functie == "rastrigine"
			x2=range(bounds_lower[1],bounds_upper[1], step=0.5)
			y2=range(bounds_lower[2],bounds_upper[2], step=0.5)
			A=10
			d=2
		    f(x2,y2) = d*A + x2^2-A*cos(2pi*x2) + y2^2-A*cos(2pi*y2)
		end
end

# ╔═╡ 9e2b4e60-42ee-11eb-0d7f-c1faa8426796
begin

	@gif for iteration in 1:T

		x_contour = []; y_contour = [];
		for bee in populations[iteration]
			append!(x_contour,bee[1]); append!(y_contour, bee[2])
		end

		plot(x2,y2,f,st=:contour,
		label="Objective function",
		xlims=(bounds_lower[1],bounds_upper[1]),
		ylims=(bounds_lower[2],bounds_upper[2]),
		legend=:outerbottom)

		scatter!(x_contour, y_contour,
		xlabel="x1",
		ylabel="x2",
		zlabel="x3",
		title="Evolution of populations over time",
		titlefont = font(15),
		c="blue",
		markershape=  :circle,
		label="Position of bees after iteration "*string(iteration),
		legend = :outerbottom)
		sleep(2)
	end every 1

end

# ╔═╡ 71321ef0-42eb-11eb-0635-b1ce95226c75
begin
	@gif for iteration in 1:T

		x_surface = []; y_surface = []; z_surface = []
		for bee in populations[iteration]
			append!(x_surface,bee[1]); append!(y_surface, bee[2]); append!(z_surface, 0)
		end

	plot(x2,y2,f,st=:surface,
		label="Objective function",
		# camera=(-30,30),
		xlims=(bounds_lower[1],bounds_upper[1]),
		ylims=(bounds_lower[2],bounds_upper[2]),
		zlims=zlims,
		legend=:outerbottom) #,c=my_cg) #,camera=(-30,30))

	scatter!(x_surface, y_surface, z_surface,
		xlabel="x1",
		ylabel="x2",
		# title="Evolution of populations over time",
		# titlefont = font(15),
		c="blue",
		markershape=  :circle,
		label="Position of bees after iteration "*string(iteration),
		legend = :outerbottom)
	sleep(2)
	end every 1
end

# ╔═╡ 65bf09be-4377-11eb-2415-3b8be310a065
md" With the chosen parameters, the **Sphere, Ackley** and **Rastrigine** function have a minimizer at (0,0).

The **Rosenbrock** function has a minimum value of zero at (0,0) and the **Branin** function has 3 global minima: (-pi, 12.275), (pi, 2.275) and (9.425, 2.475).  "

# ╔═╡ 4219c780-4381-11eb-2289-316eb02b282f
md"

#### Evolution of fitness during optimization

Below, for each iteration the **fitness** value for the best food source so far at that iteration is plotted."

# ╔═╡ 076d2e10-4381-11eb-3e12-6f9d9abe7f9a
plot(best_fitness_tracker,label="Fitness",	xlabel="iteration",
		ylabel="fitness", title="Evolution of fitness", legend=:outerbottom)

# ╔═╡ 749ca5a0-43a5-11eb-3ad6-453463f18aed

# ╔═╡ 4ef00630-43a5-11eb-2095-b3cde3592af9
md"#### Comparison with optim package
Below, the performance of the ABC algorithm is compared with the optim package."

# ╔═╡ 36031b80-445e-11eb-0c7c-bd597ac43967
f_optimize

# ╔═╡ 895b6c20-444e-11eb-0db8-3f6ce2a9db0b
md" ###### Compare with Nelder-Mead algorithm
Nelder-Mead is the standard algorithm when no gradient information is available. Instead, Nelder-Mead is a direct search method that keeps track of the function value at a number of points (forming a simplex) in the search space. Given a simplex, four actions (reflect, expand, contract, or shrink) are executed to iteratively replace the worst point with a better point.

Below, the Nelder-Mead algorithm is run for the same amount of iterations as the ABC algorithm. "

# ╔═╡ 8f7a52d0-43a7-11eb-3351-eb765609e56e
result = optimize(f_optimize, starting_solution,  method =NelderMead(), iterations = T)

# ╔═╡ 4417cb40-444e-11eb-1f2d-933616b91239
Optim.minimizer(result)

# ╔═╡ 6dea1f90-444e-11eb-1501-f99dfba633ea
Optim.minimum(result)

# ╔═╡ 198e69b0-4453-11eb-304f-43a227316038
md" We see that the Nelder-Mead algorithm of the Optim package happens to get stuck in a local mininimum when trying to minimize the Ackley and Rastrigine function, while the ABC algorithm always succeeds in finding the global minimizer.

For the other test functions, both the Nelder-Mead algorithm and the ABC algorithm succeed in finding the global minimum (or one of the global minima for the Branin test function).
 "

# ╔═╡ 0c78a8b0-4379-11eb-20dc-3dd46dc183d3
md" ## References
Karaboga, D., & Basturk, B. (2007). A powerful and efficient algorithm for numerical function optimization: artificial bee colony (ABC) algorithm. Journal of global optimization, 39(3), 459-471.

"

# ╔═╡ Cell order:
# ╟─eee523f0-436e-11eb-2e92-d59ed5c533e6
# ╟─3c1c9be0-4373-11eb-1b1e-3d3f2e6afb19
# ╟─7a2fdff0-4373-11eb-337a-8fb102ef7b78
# ╠═935e57b0-438f-11eb-0ed4-b3ad455bbd63
# ╟─16e1ad62-439c-11eb-1b3e-577f23587233
# ╟─02f4f820-439c-11eb-0c6e-7509dfc20674
# ╠═0ca837e0-42ef-11eb-17fa-9335cb9a3997
# ╟─a9462a5e-4373-11eb-3b48-39a2596229a9
# ╟─27f302ee-42ea-11eb-2d9e-49dffc0d983d
# ╟─3235a8d2-42ea-11eb-1fe1-6d91eca83dad
# ╟─085c34e0-4374-11eb-1ba7-7fb1af1d38a4
# ╠═0388b3ce-4374-11eb-03f4-9b4c74bd5ff0
# ╟─6bd962e0-439c-11eb-0ce9-c348bd45b225
# ╟─350644d0-4375-11eb-2b68-8fb8f4bd7c2a
# ╟─c7c64720-437f-11eb-15a2-477ab2fe0792
# ╟─caa5b6d0-4373-11eb-0cdd-1961e6698727
# ╠═f347e610-42a3-11eb-2116-ef50f1246cf3
# ╠═54c02380-42a4-11eb-0240-7b2d895cb337
# ╟─9bb28f50-4374-11eb-2b10-e5effcbc8438
# ╟─b81d7f30-42a5-11eb-27ce-f1cc849ffdc5
# ╟─ae391ef0-47d4-11eb-15ae-a7d5f188644b
# ╠═9e2b4e60-42ee-11eb-0d7f-c1faa8426796
# ╠═71321ef0-42eb-11eb-0635-b1ce95226c75
# ╟─65bf09be-4377-11eb-2415-3b8be310a065
# ╟─4219c780-4381-11eb-2289-316eb02b282f
# ╠═076d2e10-4381-11eb-3e12-6f9d9abe7f9a
# ╟─749ca5a0-43a5-11eb-3ad6-453463f18aed
# ╟─4ef00630-43a5-11eb-2095-b3cde3592af9
# ╠═91d40c30-43a5-11eb-1805-ef48b11c3a40
# ╠═36031b80-445e-11eb-0c7c-bd597ac43967
# ╠═df912bb0-43a5-11eb-3bd2-ebc17c073757
# ╟─895b6c20-444e-11eb-0db8-3f6ce2a9db0b
# ╠═8f7a52d0-43a7-11eb-3351-eb765609e56e
# ╠═4417cb40-444e-11eb-1f2d-933616b91239
# ╠═6dea1f90-444e-11eb-1501-f99dfba633ea
# ╟─198e69b0-4453-11eb-304f-43a227316038
# ╟─0c78a8b0-4379-11eb-20dc-3dd46dc183d3

max_stars_repo_name: kirstvh/STMOZOO

max_issues_repo_name: kirstvh/STMOZOO

max_forks_repo_name: kirstvh/STMOZOO


text:
@testset "fsm_loopback.jl" begin
    TIMEOUT_DFLT = 1000
    @testset "start in LAST_ACK, ack" begin
        conn = fsm_in_last_ack()
        tick!(conn, 4TIMEOUT_DFLT)
        expect_state(conn, JLSponge.LAST_ACK)
        send_ack!(conn, WrappingInt32(2), WrappingInt32(2))
        tick!(conn, 1)
        expect_state(conn, JLSponge.CLOSED)
    end

    @testset "start in CLOSE_WAIT, close(), throw away first FIN, ack re-tx FIN" begin
        conn = fsm_in_close_wait()
        tick!(conn, 4TIMEOUT_DFLT)
        expect_state(conn, JLSponge.CLOSE_WAIT)
        end_input_stream!(conn)
        tick!(conn, 1)
        expect_state(conn, JLSponge.LAST_ACK)

        seg1 = expect_one_seg(conn; fin=true)

        tick!(conn, TIMEOUT_DFLT - 2)
        expect_no_seg(conn)
        tick!(conn, 2)

        seg2 = expect_one_seg(conn; fin=true, seqno=seg1.header.seqno)

        rx_seqno = WrappingInt32(2)
        ack_expect = rx_seqno
        send_ack!(conn, ack_expect, seg2.header.seqno - 1)
        expect_state(conn, JLSponge.LAST_ACK)

        send_ack!(conn, ack_expect, seg2.header.seqno + 1)

        tick!(conn, 1)
        expect_state(conn, JLSponge.CLOSED)
    end

    @testset "start in ESTABLISHED, send FIN, recv ACK, check for CLOSE_WAIT" begin
        conn = fsm_in_established()
        tick!(conn, 4TIMEOUT_DFLT)
        expect_state(conn, JLSponge.ESTABLISHED)

        rx_seqno = WrappingInt32(1)
        ack_expect = rx_seqno + 1
        send_fin!(conn, rx_seqno, WrappingInt32(0))
        tick!(conn, 1)

        expect_one_seg(conn, ack=true, ackno=ack_expect)
        expect_state(conn, JLSponge.CLOSE_WAIT)

        send_fin!(conn, rx_seqno, WrappingInt32(0))
        tick!(conn, 1)

        expect_one_seg(conn; ack=true, ackno=ack_expect)

        expect_state(conn, JLSponge.CLOSE_WAIT)

        tick!(conn, 1)
        end_input_stream!(conn)
        tick!(conn, 1)
        expect_state(conn, JLSponge.LAST_ACK)
        seg3 = expect_one_seg(conn; fin=true)
        send_ack!(conn, ack_expect, seg3.header.seqno + 1)
        tick!(conn, 1)
        expect_state(conn, JLSponge.CLOSED)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using miniWeather

using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    choices_weather_types = ["collision", "thermal", "mountain_waves",
        "turbulence", "density_current", "injection"]
    @add_arg_table s begin
        "--nx_glob"
        help = "Total number of cells in the x-direction"
        arg_type = Int
        default = 100
        "--nz_glob"
        help = "Total number of cells in the z-direction"
        arg_type = Int
        default = 50
        "--sim_time"
        help = "How many seconds to run the simulation"
        arg_type = Int
        default = 50
        "--output_freq"
        help = "How frequently to output data to file (in seconds)"
        arg_type = Int
        default = 10
        "--weather_type"
        help = "Choose weather scenario and how to initialize the data. " *
                "Must be one of " * join(choices_weather_types, ", ", " or ")
        range_tester = (x->x ∈ choices_weather_types)
        default = "thermal"
    end

    return parse_args(s)
end

function main()

    config = parse_commandline()

    println("Parsed args:")
    for (arg, val) in config
        println("  $arg  =>  $val")
    end

    model, grid = init(config["nx_glob"], config["nz_glob"], config["sim_time"], config["weather_type"])

    mass0, te0 = reductions(model, grid)
    println("mass = $mass0, te = $te0")

    run!(model, grid, config["output_freq"], verbose=true)

    mass, te = reductions(model, grid)
    println("d_mass: ", (mass - mass0) / mass0)
    println("d_te:   ", (te - te0) / te0)

end

main()


text:
import DSGE: Measurement, get_setting, n_observables, n_shocks_exogenous

"""
```
measurement(m::KrusellSmith{T}, TTT::Matrix{T}, RRR::Matrix{T},
            CCC::Vector{T}) where {T<:AbstractFloat}
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::KrusellSmith{T},
                     TTT::Matrix{T}, RRR::Matrix{T},
                     CCC::Vector{T}) where {T<:AbstractFloat}
    endo      = m.endogenous_states
    exo       = m.exogenous_shocks
    obs       = m.observables

    _n_model_states = n_model_states(m)

    _n_observables = n_observables(m)
    _n_shocks_exogenous = n_shocks_exogenous(m)

    ZZ = zeros(_n_observables, _n_model_states)
    DD = zeros(_n_observables)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    # Output in log levels
    ZZ[obs[:obs_gdp], endo[:z′_t]]   .= (m[:Kstar]^m[:α])*(m[:Lstar]^(1.0-m[:α]))
    ZZ[obs[:obs_gdp], endo[:K′_t]]   .= (m[:α]/m[:Kstar])*(m[:Kstar]^m[:α])*(m[:Lstar]^(1.0-m[:α]))

    # Measurement error
    EE[obs[:obs_gdp], obs[:obs_gdp]] = m[:e_y]

    # Variance of innovations
    QQ[exo[:z_sh], exo[:z_sh]] = m[:σ_z]^2

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC .!= 0)
        DD += ZZ*((UniformScaling(1) - TTT)\CCC)
    end

    return Measurement(ZZ, DD, QQ, EE)
end


text:
# DESCRIPTION
# Closed-loop aircraft dynamics.
#
# INPUT
# - x: current state vector
# - u_curr: current control vector
#
# OUTPUT
# - xdot: system dynamics
# - out: output vector
#
# ------------------------------------------------------------------------------

function F16_DI!(x,u_curr)
    # unpack states
    VT = x[1]
    phi = x[4]
    p = x[7]
    q = x[8]
    r = x[9]
    # unpack command filter states
    p_cmd = x[NSTATES+1]
    intperr_cmd = x[NSTATES+2]
    q_cmd = x[NSTATES+3]
    intqerr_cmd = x[NSTATES+4]
    r_cmd = x[NSTATES+5]
    intrerr_cmd = x[NSTATES+6]
    δ_a = x[NSTATES+7]
    δ_e = x[NSTATES+8]
    δ_r = x[NSTATES+9]
    # control perturbations
    Δδ_a = u_curr[2]-δ_a
    Δδ_e = u_curr[3]-δ_e
    Δδ_r = u_curr[4]-δ_r
    # feedforward
    pdot_cmd = -1/τ_p*p_cmd + k_a/τ_p*Δδ_a
    qdot_cmd = -1/τ_q*q_cmd + k_e/τ_q*Δδ_e
    rdot_cmd = -1/τ_r*r_cmd + k_r/τ_r*Δδ_r
    # feedback compensation
    νp = pdot_cmd + kp_roll*(p_cmd-p) + ki_roll*(intperr_cmd)
    νq = qdot_cmd + kp_pitch*(q_cmd-q) + ki_pitch*(intqerr_cmd)
    νr = rdot_cmd + kp_yaw*(r_cmd-r) + ki_yaw*(intrerr_cmd)
    # pseudo control
    ν=[νp, νq, νr]
    # total velocity [kts]
    V=VT/1.688
    # coefficient matrix interpolation
    idx1=Int(min(max(floor(V/100)+1,1),7))
    idx2=idx1+1;
    intscale=min((V-(idx1-1)*100)/100,1.)
    CA=CATAB[:,:,idx1]+intscale*(CATAB[:,:,idx2]-CATAB[:,:,idx1])
    CBinv=CBinvTAB[:,:,idx1]+intscale*(CBinvTAB[:,:,idx2]-CBinvTAB[:,:,idx1])
    # dynamic inversion
    Δu=CBinv*(ν-CA*[p, q, r])
    # control inputs
    u=[u_curr[1]; Δu[1]+δ_a; Δu[2]+δ_e; Δu[3]+δ_r]
    # state derivative
    der, out = F16!(x[1:NSTATES],u)
    # rigid-body state derivative
    xdot=zeros(NSTATES+NDISTATES,1)
    xdot[1:NSTATES]=der
    # command filters state derivative
    xdot[NSTATES+1] = pdot_cmd
    xdot[NSTATES+2] = p_cmd-p
    xdot[NSTATES+3] = qdot_cmd+r*sin(phi)
    xdot[NSTATES+4] = q_cmd-q
    xdot[NSTATES+5] = rdot_cmd+G/VT*phi
    xdot[NSTATES+6] = r_cmd-r
    xdot[NSTATES+7] = (u_curr[2]-δ_a)/τ_u
    xdot[NSTATES+8] = (u_curr[3]-δ_e)/τ_u
    xdot[NSTATES+9] = (u_curr[4]-δ_r)/τ_u
    #
    return xdot, out
end


text:
import Base.Threads

export push_v!

function push_v!(p, kernel:: CloudInCell, m :: TwoDGrid, ex, ey, bz, dt)

    nbpart = size(p.array)[2]

    dx = m.dx
    dy = m.dy

    Threads.@threads for ipart = 1:nbpart

        v1 = p.array[3, ipart]
        v2 = p.array[4, ipart]

        xp = p.array[1, ipart] / dx
        yp = p.array[2, ipart] / dy

        i = floor(Int, xp) + 1
        j = floor(Int, yp) + 1

        dxp = xp - i + 1
        dyp = yp - j + 1

        a1 = (1 - dxp) * (1 - dyp)
        a2 = dxp * (1 - dyp)
        a3 = dxp * dyp
        a4 = (1 - dxp) * dyp

        e1 = a1 * ex[i, j] + a2 * ex[i+1, j] + a3 * ex[i+1, j+1] + a4 * ex[i, j+1]
        e2 = a1 * ey[i, j] + a2 * ey[i+1, j] + a3 * ey[i+1, j+1] + a4 * ey[i, j+1]
        b3 = a1 * bz[i, j] + a2 * bz[i+1, j] + a3 * bz[i+1, j+1] + a4 * bz[i, j+1]

        v1 += 0.5dt * e1
        v2 += 0.5dt * e2

        tantheta = 0.5dt * b3
        sintheta = 2 * tantheta / (1 + tantheta * tantheta)

        v1 += v2 * tantheta
        v2 += -v1 * sintheta
        v1 += v2 * tantheta

        p.array[3, ipart] = v1 + 0.5dt * e1
        p.array[4, ipart] = v2 + 0.5dt * e2

    end

end

export push_x!

function push_x!(p, mesh::TwoDGrid, dt::Float64)

    nbpart = size(p.array)[2]

    dimx, dimy = mesh.dimx, mesh.dimy

    Threads.@threads for i = 1:nbpart
        p1 = p.array[1, i] + dt * p.array[3, i]
        p2 = p.array[2, i] + dt * p.array[4, i]
        p1 > dimx && (p1 -= dimx)
        p2 > dimy && (p2 -= dimy)
        p1 < 0.0 && (p1 += dimx)
        p2 < 0.0 && (p2 += dimy)
        p.array[1, i] = p1
        p.array[2, i] = p2
    end

end


text:
# Autogenerated wrapper script for soapysdr_jll for aarch64-apple-darwin
export libsoapysdr, soapysdrutil

JLLWrappers.@generate_wrapper_header("soapysdr")
JLLWrappers.@declare_library_product(libsoapysdr, "@rpath/libSoapySDR.0.8.dylib")
JLLWrappers.@declare_executable_product(soapysdrutil)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libsoapysdr,
        "lib/libSoapySDR.0.8.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        soapysdrutil,
        "bin/SoapySDRUtil",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()