text:
#
# Powell.jl --
#
# Mike Powell's derivative free optimization methods for Julia.
#
# ----------------------------------------------------------------------------
#
# This file is part of OptimPack.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 2015, Éric Thiébaut.
#
# ----------------------------------------------------------------------------

# Possible status values returned by COBYLA.
const COBYLA_INITIAL_ITERATE       = convert(Cint,  2)
const COBYLA_ITERATE               = convert(Cint,  1)
const COBYLA_SUCCESS               = convert(Cint,  0)
const COBYLA_ROUNDING_ERRORS       = convert(Cint, -1)
const COBYLA_TOO_MANY_EVALUATIONS  = convert(Cint, -2)
const COBYLA_BAD_ADDRESS           = convert(Cint, -3)
const COBYLA_CORRUPTED             = convert(Cint, -4)

# Get a textual explanation of the status returned by COBYLA.
function cobyla_reason(status::Integer)
    ptr = ccall((:cobyla_reason, opklib), Ptr{UInt8}, (Cint,), status)
    if ptr == C_NULL
        error("unknown COBYLA status: ", status)
    end
    bytestring(ptr)
end

# Wrapper for the objective function in COBYLA, the actual objective
# function is provided by the client data.
function cobyla_objfun(n::Cptrdiff_t, m::Cptrdiff_t, x_::Ptr{Cdouble},
                       c_::Ptr{Cdouble}, f_::Ptr{Any})
    x = pointer_to_array(x_, n)
    f = unsafe_pointer_to_objref(f_)
    fx::Cdouble = (m > 0 ? f(x, pointer_to_array(c_, m)) : f(x))
    return fx
end
const cobyla_objfun_c = cfunction(cobyla_objfun, Cdouble,
                                  (Cptrdiff_t, Cptrdiff_t, Ptr{Cdouble},
                                   Ptr{Cdouble}, Ptr{Any}))

function cobyla_check(n::Integer, m::Integer, rhobeg::Real, rhoend::Real)
    if n < 2
        "bad number of variables"
    elseif m < 0
        "bad number of constraints"
    elseif rhoend < 0 || rhoend > rhobeg
        "bad trust region radius settings"
    end
end

function cobyla!(f::Function, x::Vector{Cdouble},
                 m::Integer, rhobeg::Real, rhoend::Real;
                 verbose::Integer=0, maxeval::Integer=500)
    n = length(x)
    reason = cobyla_check(n, m, rhobeg, rhoend)
    reason == nothing || error(reason)
    w = Array(Cdouble, n*(3*n + 2*m + 11) + 4*m + 6)
    iact = Array(Cptrdiff_t, m + 1)
    maxfun = Cptrdiff_t[maxeval]
    status = ccall((:cobyla, opklib), Cint,
                   (Cptrdiff_t, Cptrdiff_t, Ptr{Void}, Ptr{Any},
                    Ptr{Cdouble}, Cdouble, Cdouble, Cptrdiff_t,
                    Ptr{Cptrdiff_t}, Ptr{Cdouble}, Ptr{Cptrdiff_t}),
                   n, m, cobyla_objfun_c, pointer_from_objref(f),
                   x, rhobeg, rhoend, verbose, maxfun, w, iact)
    if status != COBYLA_SUCCESS
        error(cobyla_reason(status))
    end
    #return w[1]
end

# Context for reverse communication variant of COBYLA.
type CobylaContext
    ptr::Ptr{Void}
    n::Int
    m::Int
    rhobeg::Cdouble
    rhoend::Cdouble
    verbose::Int
    maxeval::Int
end

# Create a new reverse communication workspace for COBYLA algorithm.
# A typical usage is:
# ```
# x = Array(Cdouble, n)
# c = Array(Cdouble, m)
# x[...] = ... # initial solution
# ctx = cobyla_create(n, m, rhobeg, rhoend, verbose=1, maxeval=500)
# status = cobyla_get_status(ctx)
# while status == COBYLA_ITERATE
#     fx = ...       # compute function value at X
#     c[...] = ...   # compute constraints at X
#     status = cobyla_iterate(ctx, fx, x, c)
# end
# if status != COBYLA_SUCCESS
#     println("Something wrong occured in COBYLA: ", cobyla_reason(status))
# end
# ```
function cobyla_create(n::Integer, m::Integer,
                       rhobeg::Real, rhoend::Real;
                       verbose::Integer=0, maxeval::Integer=500)
    reason = cobyla_check(n, m, rhobeg, rhoend)
    reason == nothing || error(reason)
    ptr = ccall((:cobyla_create, opklib), Ptr{Void},
                (Cptrdiff_t, Cptrdiff_t, Cdouble, Cdouble,
                 Cptrdiff_t, Cptrdiff_t),
                n, m, rhobeg, rhoend, verbose, maxeval)
    if ptr == C_NULL
        reason = (errno() == Base.Errno.ENOMEM
                  ? "insufficient memory"
                  : "unexpected error")
        error(reason)
    end
    ctx = CobylaContext(ptr, n, m, rhobeg, rhoend, verbose, maxeval)
    finalizer(ctx, ctx -> ccall((:cobyla_delete, opklib), Void,
                                (Ptr{Void},), ctx.ptr))
    return ctx
end

# Perform the next iteration of the reverse communication variant of the
# COBYLA algorithm.  On entry, the workspace status must be
# `COBYLA_ITERATE`, `f` and `c` are the function value and the constraints
# at `x`.  On exit, the returned value (the new workspace status) is:
# `COBYLA_ITERATE` if a new trial point has been stored in `x` and if user
# is requested to compute the function value and the constraints on the new
# point; `COBYLA_SUCCESS` if algorithm has converged and `x` has been set
# with the variables at the solution (the corresponding function value can
# be retrieved with `cobyla_get_last_f`); anything else indicates an error
# (see `cobyla_reason` for an explanatory message).
function cobyla_iterate(ctx::CobylaContext, f::Real, x::Vector{Cdouble},
                        c::Vector{Cdouble})
    length(x) == ctx.n || error("bad number of variables")
    length(c) == ctx.m || error("bad number of constraints")
    ccall((:cobyla_iterate, opklib), Cint,
          (Ptr{Void}, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
          ctx.ptr, f, x, c)
end
# The same but without constraints.
function cobyla_iterate(ctx::CobylaContext, f::Real, x::Vector{Cdouble})
    length(x) == ctx.n || error("bad number of variables")
    ctx.m == 0 || error("bad number of constraints")
    ccall((:cobyla_iterate, opklib), Cint,
          (Ptr{Void}, Cdouble, Ptr{Cdouble}, Ptr{Void}),
          ctx.ptr, f, x, C_NULL)
end

# Restart COBYLA algorithm using the same parameters.  The return value is
# the new status of the algorithm, see `cobyla_get_status` for details.
function cobyla_restart(ctx::CobylaContext)
    ccall((:cobyla_restart, opklib), Cint, (Ptr{Void},), ctx.ptr)
end

# Get the current status of the algorithm.  Result is: `COBYLA_ITERATE` if
# user is requested to compute F(X) and C(X); `COBYLA_SUCCESS` if algorithm
# has converged; anything else indicates an error (see `cobyla_reason` for
# an explanatory message).
function cobyla_get_status(ctx::CobylaContext)
    ccall((:cobyla_get_status, opklib), Cint, (Ptr{Void},), ctx.ptr)
end

# Get the current number of function evaluations.  Result is -1 if
# something is wrong (e.g. CTX is NULL), nonnegative otherwise.
function cobyla_get_nevals(ctx::CobylaContext)
    ccall((:cobyla_get_nevals, opklib), Cptrdiff_t, (Ptr{Void},), ctx.ptr)
end

# Get the current size of the trust region.  Result is 0 if algorithm not
# yet started (before first iteration), -1 if something is wrong (e.g. CTX
# is NULL), strictly positive otherwise.
function cobyla_get_rho(ctx::CobylaContext)
    ccall((:cobyla_get_rho, opklib), Cdouble, (Ptr{Void},), ctx.ptr)
end

# Get the last function value.  Upon convergence of `cobyla_iterate`
# (i.e. return with status `COBYLA_SUCCESS`), this value corresponds to the
# function at the solution; otherwise, this value corresponds to the
# previous set of variables.
function cobyla_get_last_f(ctx::CobylaContext)
    ccall((:cobyla_get_last_f, opklib), Cdouble, (Ptr{Void},), ctx.ptr)
end

# Get a textual explanation of the current status.
cobyla_get_reason(ctx::CobylaContext) = cobyla_reason(cobyla_get_status(ctx))

function cobyla_test(revcom::Bool=false)
    # Beware that order of operations may affect the result (whithin
    # rounding errors).  I have tried to keep the same ordering as F2C
    # which takes care of that, in particular when converting expressions
    # involving powers.
    prt(s) = println("\n       "*s)
    for nprob in 1:10
        if nprob == 1
            # Minimization of a simple quadratic function of two variables.
            prt("Output from test problem 1 (Simple quadratic)")
            n = 2
            m = 0
            xopt = Array(Cdouble, n)
            xopt[1] = -1.0
            xopt[2] = 0.0
            function ftest(x::Vector{Cdouble})
                r1 = x[1] + 1.0
                r2 = x[2]
                fc = 10.0*(r1*r1) + (r2*r2)
                return fc
            end
        elseif nprob == 2
            # Easy two dimensional minimization in unit circle.
            prt("Output from test problem 2 (2D unit circle calculation)")
            n = 2
            m = 1
            xopt = Array(Cdouble, n)
            xopt[1] = sqrt(0.5)
            xopt[2] = -xopt[1]
            function ftest(x::Vector{Cdouble}, con::Vector{Cdouble})
                fc = x[1]*x[2]
                con[1] = 1.0 - x[1]*x[1] - x[2]*x[2]
                return fc
            end
        elseif nprob == 3
            # Easy three dimensional minimization in ellipsoid.
            prt("Output from test problem 3 (3D ellipsoid calculation)")
            n = 3
            m = 1
            xopt = Array(Cdouble, n)
            xopt[1] = 1.0/sqrt(3.0)
            xopt[2] = 1.0/sqrt(6.0)
            xopt[3] = -0.33333333333333331
            function ftest(x::Vector{Cdouble}, con::Vector{Cdouble})
                fc = x[1]*x[2]*x[3]
                con[1] = 1.0 - (x[1]*x[1]) - 2.0*(x[2]*x[2]) - 3.0*(x[3]*x[3])
                return fc
            end
        elseif nprob == 4
            # Weak version of Rosenbrock's problem.
            prt("Output from test problem 4 (Weak Rosenbrock)")
            n = 2
            m = 0
            xopt = Array(Cdouble, n)
            xopt[1] = -1.0
            xopt[2] = 1.0
            function ftest(x::Vector{Cdouble})
                r2 = x[1]
                r1 = r2*r2 - x[2]
                r3 = x[1] + 1.0
                fc = r1*r1 + r3*r3
                return fc
            end
        elseif nprob == 5
            # Intermediate version of Rosenbrock's problem.
            prt("Output from test problem 5 (Intermediate Rosenbrock)")
            n = 2
            m = 0
            xopt = Array(Cdouble, n)
            xopt[1] = -1.0
            xopt[2] = 1.0
            function ftest(x::Vector{Cdouble})
                r2 = x[1]
                r1 = r2*r2 - x[2]
                r3 = x[1] + 1.0
                fc = r1*r1*10.0 + r3*r3
                return fc
            end
        elseif nprob == 6
            # This problem is taken from Fletcher's book Practical Methods
            # of Optimization and has the equation number (9.1.15).
            prt("Output from test problem 6 (Equation (9.1.15) in Fletcher)")
            n = 2
            m = 2
            xopt = Array(Cdouble, n)
            xopt[1] = sqrt(0.5)
            xopt[2] = xopt[1]
            function ftest(x::Vector{Cdouble}, con::Vector{Cdouble})
                fc = -x[1] - x[2]
                r1 = x[1]
                con[1] = x[2] - r1*r1
                r1 = x[1]
                r2 = x[2]
                con[2] = 1.0 - r1*r1 - r2*r2
                return fc
            end
        elseif nprob == 7
            # This problem is taken from Fletcher's book Practical Methods
            # of Optimization and has the equation number (14.4.2).
            prt("Output from test problem 7 (Equation (14.4.2) in Fletcher)")
            n = 3
            m = 3
            xopt = Array(Cdouble, n)
            xopt[1] = 0.0
            xopt[2] = -3.0
            xopt[3] = -3.0
            function ftest(x::Vector{Cdouble}, con::Vector{Cdouble})
                fc = x[3]
                con[1] = x[1]*5.0 - x[2] + x[3]
                r1 = x[1]
                r2 = x[2]
                con[2] = x[3] - r1*r1 - r2*r2 - x[2]*4.0
                con[3] = x[3] - x[1]*5.0 - x[2]
                return fc
            end
        elseif nprob == 8
            # This problem is taken from page 66 of Hock and Schittkowski's
            # book Test Examples for Nonlinear Programming Codes. It is
            # their test problem Number 43, and has the name Rosen-Suzuki.
            prt("Output from test problem 8 (Rosen-Suzuki)")
            n = 4
            m = 3
            xopt = Array(Cdouble, n)
            xopt[1] = 0.0
            xopt[2] = 1.0
            xopt[3] = 2.0
            xopt[4] = -1.0
            function ftest(x::Vector{Cdouble}, con::Vector{Cdouble})
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                r4 = x[4]
                fc = (r1*r1 + r2*r2 + r3*r3*2.0 + r4*r4 - x[1]*5.0
                      - x[2]*5.0 - x[3]*21.0 + x[4]*7.0)
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                r4 = x[4]
                con[1] = (8.0 - r1*r1 - r2*r2 - r3*r3 - r4*r4 - x[1]
                          + x[2] - x[3] + x[4])
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                r4 = x[4]
                con[2] = (10.0 - r1*r1 - r2*r2*2.0 - r3*r3 - r4*r4*2.0
                          + x[1] + x[4])
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                con[3] = (5.0 - r1*r1*2.0 - r2*r2 - r3*r3 - x[1]*2.0
                          + x[2] + x[4])
                return fc
            end
        elseif nprob == 9
            # This problem is taken from page 111 of Hock and
            # Schittkowski's book Test Examples for Nonlinear Programming
            # Codes. It is their test problem Number 100.
            prt("Output from test problem 9 (Hock and Schittkowski 100)")
            n = 7
            m = 4
            xopt = Array(Cdouble, n)
            xopt[1] =  2.330499
            xopt[2] =  1.951372
            xopt[3] = -0.4775414
            xopt[4] =  4.365726
            xopt[5] = -0.624487
            xopt[6] =  1.038131
            xopt[7] =  1.594227
            function ftest(x::Vector{Cdouble}, con::Vector{Cdouble})
                r1 = x[1] - 10.0
                r2 = x[2] - 12.0
                r3 = x[3]
                r3 *= r3
                r4 = x[4] - 11.0
                r5 = x[5]
                r5 *= r5
                r6 = x[6]
                r7 = x[7]
                r7 *= r7
                fc = (r1*r1 + r2*r2*5.0 + r3*r3 + r4*r4*3.0
                      + r5*(r5*r5)*10.0 + r6*r6*7.0 + r7*r7
                      - x[6]*4.0*x[7] - x[6]*10.0 - x[7]*8.0)
                r1 = x[1]
                r2 = x[2]
                r2 *= r2
                r3 = x[4]
                con[1] = (127.0 - r1*r1*2.0 - r2*r2*3.0 - x[3]
                          - r3*r3*4.0 - x[5]*5.0)
                r1 = x[3]
                con[2] = (282.0 - x[1]*7.0 - x[2]*3.0 - r1*r1*10.0
                          - x[4] + x[5])
                r1 = x[2]
                r2 = x[6]
                con[3] = (196.0 - x[1]*23.0 - r1*r1 - r2*r2*6.0
                          + x[7]*8.0)
                r1 = x[1]
                r2 = x[2]
                r3 = x[3]
                con[4] = (r1*r1*-4.0 - r2*r2 + x[1]*3.0*x[2]
                          - r3*r3*2.0 - x[6]*5.0 + x[7]*11.0)
                return fc
            end
        elseif nprob == 10
            # This problem is taken from page 415 of Luenberger's book
            # Applied Nonlinear Programming. It is to maximize the area of
            # a hexagon of unit diameter.
            prt("Output from test problem 10 (Hexagon area)")
            n = 9
            m = 14
            xopt = Array(Cdouble, n)
            function ftest(x::Vector{Cdouble}, con::Vector{Cdouble})
                fc = -0.5*(x[1]*x[4] - x[2]*x[3] + x[3]*x[9] - x[5]*x[9]
                           + x[5]*x[8] - x[6]*x[7])
                r1 = x[3]
                r2 = x[4]
                con[1] = 1.0 - r1*r1 - r2*r2
                r1 = x[9]
                con[2] = 1.0 - r1*r1
                r1 = x[5]
                r2 = x[6]
                con[3] = 1.0 - r1*r1 - r2*r2
                r1 = x[1]
                r2 = x[2] - x[9]
                con[4] = 1.0 - r1*r1 - r2*r2
                r1 = x[1] - x[5]
                r2 = x[2] - x[6]
                con[5] = 1.0 - r1*r1 - r2*r2
                r1 = x[1] - x[7]
                r2 = x[2] - x[8]
                con[6] = 1.0 - r1*r1 - r2*r2
                r1 = x[3] - x[5]
                r2 = x[4] - x[6]
                con[7] = 1.0 - r1*r1 - r2*r2
                r1 = x[3] - x[7]
                r2 = x[4] - x[8]
                con[8] = 1.0 - r1*r1 - r2*r2
                r1 = x[7]
                r2 = x[8] - x[9]
                con[9] = 1.0 - r1*r1 - r2*r2
                con[10] = x[1]*x[4] - x[2]*x[3]
                con[11] = x[3]*x[9]
                con[12] = -x[5]*x[9]
                con[13] = x[5]*x[8] - x[6]*x[7]
                con[14] = x[9]
                return fc
            end
        else
            error("bad problem number ($nprob)")
        end

        x = Array(Cdouble, n)
        for icase in 1:2
            for i in 1:n
                x[i] = 1.0
            end
            rhobeg = 0.5
            rhoend = (icase == 2 ? 1e-4 : 0.001)
            if revcom
                # Test the reverse communication variant.
                c = (m > 0 ? Array(Cdouble, m) : nothing)
                ctx = cobyla_create(n, m, rhobeg, rhoend, verbose=1, maxeval=2000)
                status = cobyla_get_status(ctx)
                while status == COBYLA_ITERATE
                    if m > 0
                        # Some constraints.
                        fx = ftest(x, c)
                        status = cobyla_iterate(ctx, fx, x, c)
                    else
                        # No constraints.
                        fx = ftest(x)
                        status = cobyla_iterate(ctx, fx, x)
                    end
                end
                if status != COBYLA_SUCCESS
                    println("Something wrong occured in COBYLA: ", cobyla_reason(status))
                end
            else
                cobyla!(ftest, x, m, rhobeg, rhoend, verbose=1, maxeval=2000)
            end
            if nprob == 10
                tempa = x[1] + x[3] + x[5] + x[7]
                tempb = x[2] + x[4] + x[6] + x[8]
                tempc = 0.5/sqrt(tempa*tempa + tempb*tempb)
                tempd = tempc*sqrt(3.0)
                xopt[1] = tempd*tempa + tempc*tempb
                xopt[2] = tempd*tempb - tempc*tempa
                xopt[3] = tempd*tempa - tempc*tempb
                xopt[4] = tempd*tempb + tempc*tempa
                for i in 1:4
                    xopt[i + 4] = xopt[i]
                end
            end
            temp = 0.0
            for i in 1:n
                r1 = x[i] - xopt[i]
                temp += r1*r1
            end
            @printf("\n     Least squares error in variables =%16.6E\n", sqrt(temp))
        end
        @printf("  ------------------------------------------------------------------\n")
    end
end

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

# Possible status values returned by NEWUOA.
const NEWUOA_INITIAL_ITERATE      = convert(Cint,  2)
const NEWUOA_ITERATE              = convert(Cint,  1)
const NEWUOA_SUCCESS              = convert(Cint,  0)
const NEWUOA_BAD_NPT              = convert(Cint, -1)
const NEWUOA_ROUNDING_ERRORS      = convert(Cint, -2)
const NEWUOA_TOO_MANY_EVALUATIONS = convert(Cint, -3)
const NEWUOA_STEP_FAILED          = convert(Cint, -4)
const NEWUOA_BAD_ADDRESS          = convert(Cint, -5)
const NEWUOA_CORRUPTED            = convert(Cint, -6)

# Get a textual explanation of the status returned by NEWUOA.
function newuoa_reason(status::Integer)
    ptr = ccall((:newuoa_reason, opklib), Ptr{UInt8}, (Cint,), status)
    if ptr == C_NULL
        error("unknown NEWUOA status: ", status)
    end
    bytestring(ptr)
end

# Wrapper for the objective function in NEWUOA, the actual objective
# function is provided by the client data.
function newuoa_objfun(n::Cptrdiff_t, x_::Ptr{Cdouble}, f_::Ptr{Any})
    x = pointer_to_array(x_, n)
    f = unsafe_pointer_to_objref(f_)
    fx::Cdouble = f(x)
    return fx
end
const newuoa_objfun_c = cfunction(newuoa_objfun, Cdouble,
                                  (Cptrdiff_t, Ptr{Cdouble}, Ptr{Any}))

function newuoa_check(n::Integer, npt::Integer, rhobeg::Real, rhoend::Real)
    if n < 2
        "bad number of variables"
    elseif npt < n + 2 || npt > div((n + 2)*(n + 1),2)
        "NPT is not in the required interval"
    elseif rhoend < 0 || rhoend > rhobeg
        "bad trust region radius settings"
    end
end

function newuoa!(f::Function, x::Vector{Cdouble},
                 npt::Integer, rhobeg::Real, rhoend::Real;
                 verbose::Integer=0, maxeval::Integer=500)
    n = length(x)
    reason = newuoa_check(n, npt, rhobeg, rhoend)
    reason == nothing || error(reason)
    w = Array(Cdouble, (npt + 13)*(npt + n) + div(3*n*(n + 3),2))
    status = ccall((:newuoa, opklib), Cint,
                   (Cptrdiff_t, Cptrdiff_t, Ptr{Void}, Ptr{Any},
                    Ptr{Cdouble}, Cdouble, Cdouble, Cptrdiff_t,
                    Cptrdiff_t, Ptr{Cdouble}),
                   n, npt, newuoa_objfun_c, pointer_from_objref(f),
                   x, rhobeg, rhoend, verbose, maxeval, w)
    if status != NEWUOA_SUCCESS
        error(newuoa_reason(status))
    end
    return nothing
end

# Context for reverse communication variant of NEWUOA.
type NewuoaContext
    ptr::Ptr{Void}
    n::Int
    npt::Int
    rhobeg::Cdouble
    rhoend::Cdouble
    verbose::Int
    maxeval::Int
end

# Create a new reverse communication workspace for NEWUOA algorithm.
# A typical usage is:
# ```
# x = Array(Cdouble, n)
# x[...] = ... # initial solution
# ctx = newuoa_create(n, npt, rhobeg, rhoend, verbose=0, maxeval=500)
# status = newuoa_get_status(ctx)
# while status == NEWUOA_ITERATE
#   fx = ... # compute function value at X
#   status = newuoa_iterate(ctx, fx, x)
# end
# if status != NEWUOA_SUCCESS
#   println("Something wrong occured in NEWUOA: ", newuoa_reason(status))
# end
# ```
function newuoa_create(n::Integer, npt::Integer,
                       rhobeg::Real, rhoend::Real;
                       verbose::Integer=0, maxeval::Integer=500)
    reason = newuoa_check(n, npt, rhobeg, rhoend)
    reason == nothing || error(reason)
    ptr = ccall((:newuoa_create, opklib), Ptr{Void},
                (Cptrdiff_t, Cptrdiff_t, Cdouble, Cdouble,
                 Cptrdiff_t, Cptrdiff_t),
                n, npt, rhobeg, rhoend, verbose, maxeval)
    if ptr == C_NULL
        reason = (errno() == Base.Errno.ENOMEM
                  ? "insufficient memory"
                  : "unexpected error")
        error(reason)
    end
    ctx = NewuoaContext(ptr, n, npt, rhobeg, rhoend, verbose, maxeval)
    finalizer(ctx, ctx -> ccall((:newuoa_delete, opklib), Void,
                                (Ptr{Void},), ctx.ptr))
    return ctx
end

# Perform the next iteration of the reverse communication version of the
# NEWUOA algorithm.  On entry, the wokspace status must be
# `NEWUOA_ITERATE`, `f` is the function value at `x`.  On exit, the
# returned value (the new wokspace status) is: `NEWUOA_ITERATE` if a new
# trial point has been stored in `x` and if user is requested to compute
# the function value for the new point; `NEWUOA_SUCCESS` if algorithm has
# converged; anything else indicates an error (see `newuoa_reason` for an
# explanatory message).
function newuoa_iterate(ctx::NewuoaContext, f::Real, x::Vector{Cdouble})
    length(x) == ctx.n || error("bad number of variables")
    ccall((:newuoa_iterate, opklib), Cint,
          (Ptr{Void}, Cdouble, Ptr{Cdouble}),
          ctx.ptr, f, x)
end

# Restart NEWUOA algorithm using the same parameters.  The return value is the
# new status of the algorithm, see `newuoa_get_status` for details.
function newuoa_restart(ctx::NewuoaContext)
    ccall((:newuoa_restart, opklib), Cint, (Ptr{Void},), ctx.ptr)
end

# Get the current status of the algorithm.  Result is: `NEWUOA_ITERATE` if
# user is requested to compute F(X); `NEWUOA_SUCCESS` if algorithm has
# converged; anything else indicates an error (see `newuoa_reason` for an
# explanatory message).
function newuoa_get_status(ctx::NewuoaContext)
    ccall((:newuoa_get_status, opklib), Cint, (Ptr{Void},), ctx.ptr)
end

# Get the current number of function evaluations.  Result is -1 if
# something is wrong (e.g. CTX is NULL), nonnegative otherwise.
function newuoa_get_nevals(ctx::NewuoaContext)
    ccall((:newuoa_get_nevals, opklib), Cptrdiff_t, (Ptr{Void},), ctx.ptr)
end

# Get the current size of the trust region.  Result is 0 if algorithm not
# yet started (before first iteration), -1 if something is wrong (e.g. CTX
# is NULL), strictly positive otherwise.
function newuoa_get_rho(ctx::NewuoaContext)
    ccall((:newuoa_get_rho, opklib), Cdouble, (Ptr{Void},), ctx.ptr)
end

function newuoa_test(revcom::Bool=false)
    # The Chebyquad test problem (Fletcher, 1965) for N = 2,4,6 and 8, with
    # NPT = 2N+1.
    function ftest(x::Vector{Cdouble})
        n = length(x)
        np = n + 1
        y = Array(Cdouble, np, n)
        for j in 1:n
            y[1,j] = 1.0
            y[2,j] = x[j]*2.0 - 1.0
        end
        for i in 2:n
            for j in 1:n
                y[i+1,j] = y[2,j]*2.0*y[i,j] - y[i-1,j]
            end
        end
        f = 0.0
        iw = 1
        for i in 1:np
            sum = 0.0
            for j in 1:n
                sum += y[i,j]
            end
            sum /= n
            if iw > 0
                sum += 1.0/(i*i - 2*i)
            end
            iw = -iw
            f += sum*sum
        end
        return f
    end

    # Run the tests.
    rhoend = 1e-6
    for n = 2:2:8
        npt = 2*n + 1
        x = Array(Cdouble, n)
        for i in 1:n
            x[i] = i/(n + 1)
        end
        rhobeg = x[1]*0.2
        @printf("\n\n    Results with N =%2d and NPT =%3d\n", n, npt)
        if revcom
            # Test the reverse communication variant.
            ctx = newuoa_create(n, npt, rhobeg, rhoend, verbose=2, maxeval=5000)
            status = newuoa_get_status(ctx)
            while status == NEWUOA_ITERATE
                fx = ftest(x)
                status = newuoa_iterate(ctx, fx, x)
            end
            if status != NEWUOA_SUCCESS
                println("Something wrong occured in NEWUOA: ", newuoa_reason(status))
            end
        else
            newuoa!(ftest, x, npt, rhobeg, rhoend, verbose=2, maxeval=5000)
        end
    end
end

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

const BOBYQA_SUCCESS              = convert(Cint,  0)
const BOBYQA_BAD_NPT              = convert(Cint, -1)
const BOBYQA_TOO_CLOSE            = convert(Cint, -2)
const BOBYQA_ROUNDING_ERRORS      = convert(Cint, -3)
const BOBYQA_TOO_MANY_EVALUATIONS = convert(Cint, -4)
const BOBYQA_STEP_FAILED          = convert(Cint, -5)

function bobyqa_reason(status::Integer)
    if status == BOBYQA_SUCCESS
        return "algorithm converged"
    elseif status == BOBYQA_BAD_NPT
        return "NPT is not in the required interval"
    elseif status == BOBYQA_TOO_CLOSE
        return "insufficient space between the bounds"
    elseif status == BOBYQA_ROUNDING_ERRORS
        return "too much cancellation in a denominator"
    elseif status == BOBYQA_TOO_MANY_EVALUATIONS
        return "maximum number of function evaluations exceeded"
    elseif status == BOBYQA_STEP_FAILED
        return "a trust region step has failed to reduce quadratic approximation"
    else
        return "unknown BOBYQA status"
    end
end

# Wrapper for the objective function in BOBYQA, the actual objective
# function is provided by the client data.
function bobyqa_objfun(n::Cptrdiff_t, x_::Ptr{Cdouble}, f_::Ptr{Any})
    x = pointer_to_array(x_, n)
    f = unsafe_pointer_to_objref(f_)
    fx::Cdouble = f(x)
    return fx
end
const bobyqa_objfun_c = cfunction(bobyqa_objfun, Cdouble,
                                  (Cptrdiff_t, Ptr{Cdouble}, Ptr{Any}))

function bobyqa!(f::Function, x::Vector{Cdouble},
                 xl::Vector{Cdouble}, xu::Vector{Cdouble},
                 rhobeg::Real, rhoend::Real;
                 npt::Union{Integer,Void}=nothing,
                 verbose::Integer=0, maxeval::Integer=500)
    n = length(x)
    length(xl) == n || error("bad length for inferior bound")
    length(xu) == n || error("bad length for superior bound")
    if npt == nothing
        npt = 2*n + 1
    end
    w = Array(Cdouble, (npt + 5)*(npt + n) + div(3*n*(n + 5),2))
    status = ccall((:bobyqa, opklib), Cint,
                (Cptrdiff_t, Cptrdiff_t, Ptr{Void}, Ptr{Any},
                 Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                 Cdouble, Cdouble, Cptrdiff_t, Cptrdiff_t,
                 Ptr{Cdouble}),
                   n, npt, bobyqa_objfun_c, pointer_from_objref(f),
                   x, xl, xu, rhobeg, rhoend, verbose, maxeval, w)
    if status != BOBYQA_SUCCESS
        error(bobyqa_reason(status))
    end
    return w[1]
end

function bobyqa_test()
    # The test function.
    function f(x::Vector{Cdouble})
        fx = 0.0
        n = length(x)
        for i in 4:2:n
            for j in 2:2:i-2
                tempa = x[i - 1] - x[j - 1]
                tempb = x[i] - x[j]
                temp = tempa*tempa + tempb*tempb
                temp = max(temp,1e-6)
                fx += 1.0/sqrt(temp)
            end
        end
        return fx
    end

    # Run the tests.
    bdl = -1.0
    bdu =  1.0
    rhobeg = 0.1
    rhoend = 1e-6
    for m in (5,10)
        q = 2.0*pi/m
        n = 2*m
        x = Array(Cdouble, n)
        xl = Array(Cdouble, n)
        xu = Array(Cdouble, n)
        for i in 1:n
            xl[i] = bdl
            xu[i] = bdu
        end
        for jcase in 1:2
            if jcase == 2
                npt = 2*n + 1
            else
                npt = n + 6
            end
            @printf("\n\n     2D output with M =%4ld,  N =%4ld  and  NPT =%4ld\n",
                    m, n, npt)
            for j in 1:m
                temp = q*j
                x[2*j - 1] = cos(temp)
                x[2*j]     = sin(temp)
            end
            fx = bobyqa!(f, x, xl, xu, rhobeg, rhoend, npt=npt, verbose=2, maxeval=500000)
            @printf("\n***** least function value: %.15e\n", fx)
        end
    end
end


text:
import Conda

doit = get(ENV, "INSTALL_GMT", "")

if (doit != "")

	depfile = joinpath(dirname(@__FILE__), "deps.jl")

	if Sys.iswindows()
		if !isfile("C:\\programs\\gmt6\\bin\\gmt.exe")		# If we have none, download and install from installer
			fn = download("http://fct-gmt.ualg.pt/gmt/data/wininstallers/gmt-win64.exe", "GMTinstaller.exe")
			run(`cmd /k GMTinstaller.exe /S`)
			rm(fn, force=true)
		end
		GMT_bindir = "C:\\programs\\gmt6\\bin"
	else
		Conda.add_channel("conda-forge")
		Conda.add("gmt")
		GMT_bindir = joinpath(Conda.ROOTENV, "bin")
	end

	libgmt = abspath(chop(read(`$(joinpath("$(GMT_bindir)", "gmt")) --show-library`, String)))

	@static Sys.iswindows() ? libgdal = "gdal_w64.dll" : (
		Sys.isapple() ? (libgdal = joinpath(Conda.ROOTENV, "lib", string(split(readlines(pipeline(`otool -L $(libgmt)`, `grep libgdal`))[1])[1])[8:end]) )  : (
				Sys.isunix() ? (libgdal = string(split(readlines(pipeline(`ldd $(libgmt)`, `grep libgdal`))[1])[3])) :
				error("Don't know how to install this package in this OS.")
			)
		)

	@static Sys.iswindows() ? libproj = "proj_w64.dll" : (
		Sys.isapple() ? (libproj = joinpath(Conda.ROOTENV, "lib", string(split(readlines(pipeline(`otool -L $(libgdal)`, `grep libproj`))[1])[1])[8:end]) )  : (
				Sys.isunix() ? (libproj = string(split(readlines(pipeline(`ldd $(libgdal)`, `grep libproj`))[1])[3])) :
				error("Don't know how to use PROJ4 in this OS.")
			)
		)

	# Save shared names in file so that GMT.jl can read them at pre-compile time
	open(depfile, "w") do f
		println(f, "_GMT_bindir = \"", escape_string(GMT_bindir), '"')
		println(f, "_libgmt  = \"", escape_string(libgmt), '"')
		println(f, "_libgdal = \"", escape_string(joinpath(GMT_bindir, libgdal)), '"')
		println(f, "_libproj = \"", escape_string(joinpath(GMT_bindir, libproj)), '"')
	end
end

ext: jl


text:
#=
exp_mult.jl
translated from matlab/C file mri_exp_mult_mex.c
originally by Sangwoo Lee, University of Michigan, Jun, 2004
based on the 2003 file int_tim_seg.m by Brad Sutton
2004-6-21 modified by JF
2020-02 by Daniel Wan
=#

export exp_mult

"""
    D = exp_mult(A, u, v ; warnboth)

Memory efficient and fast implementation of `D = A' * exp(-u * v^T)`
that is useful for B0-field-corrected MRI image reconstruction.

in:
* `A [N L]` matrix
* `u [N]` vector
* `v [M]` vector
* `warnboth` warn if both `u` and `v` are complex; default: true

out:
* `D [L M]` complex vector: `D = A' * exp(-u * v^T)`
`D_lm = sum_n A_nl^* exp(-u_n v_m)`
"""
function exp_mult(A, u::AbstractVector{<:Number}, v::AbstractVector{<:Number}
        ; warnboth::Bool = true)

    N, L = size(A)
    n = size(u,1)
    M = size(v,1)

    ~isreal(u) && ~isreal(v) && warnboth &&
        @warn "only one of u and v should be complex"

    T = promote_type(eltype(u), eltype(v), ComplexF32)
    D = zeros(T, L, M)

    for m = 1:M
        col = exp.(-u * v[m]) # [N] mth column of N × M matrix B = exp(-u * v^T)
        D[:,m] = A' * col
    end
    return D
end


text:
using AbstractPlotting, Observables
using StatsBase
using StatsMakie
import Colors
export trajectory_highlighter

const DEFAULT_α = 0.01

"""
    trajectory_highlighter(datasets, vals; kwargs...)
Open an interactive application for highlighting specific datasets
and properties of these datasets. `datasets` is a vector of _anything_ from
**DynamicalSystems.jl** that supports `plot_dataset` (currently `Dataset`
or `Matrix`). Each dataset corresponds to a specific value from `vals`
(a `Vector{<:Real}`). The value of `vals` gives each dataset
a specific color based on a colormap.

The application is composed of two scenes: the left scene plots the datasets,
while the right scene plots the histogram of the `vals`. The function
returns the two scenes `data_scene, hist_scene`.

## Interaction
Clicking on a bin of the histogram plot will "highlight" all data
whose value belongs in that bin. Here highlighting actually means "hidding"
(i.e. reducing their alpha value) all other data besides the ones you want
to highlight. Clicking on empty space on the histogram plot will reset
highlighting.

Clicking on a plotted series in the left window will highlight this series
as well as the histogram bin that contains its value. Clicking on empty
space will reset the highlighting.

## Keyword Arguments
* `nbins = 10, closed = :left` : used in producing the histogram.
* `α = 0.05` : the alpha value of the hidden data.
* `hα = 0.2` : the alpha value of the hidden histogram bins.
* `cmap = :viridis` : the colormap used.
* `hname = "value"` : name for the histogram axis.
* `kwargs...` : Anything else is propagated to `plot_dataset`.
"""
function trajectory_highlighter(datasets, vals;
    nbins = 10, closed=:left, α = 0.05,
    cmap = :viridis, hα = 0.2, hname = "value", kwargs...)

    N = length(datasets)
    N == length(vals) || error("data and value must have equal length")

    # First prepare the colors of the datasets:
    colormap = to_colormap(cmap, length(datasets))
    get_color(i) = Colors.color(AbstractPlotting.interpolated_getindex(
        colormap, vals[i], extrema(vals)
    ))
    # The colors are observables; the transparency can be changed
    data_αs = [Observable(1.0) for i in 1:N]
    colors = [lift(α -> RGBAf0(get_color(i), α), data_αs[i]) for i ∈ 1:N]
    data_scene = plot_datasets(datasets, colors; kwargs...)

    # now time for the histogram:
    hist = fit(StatsBase.Histogram, vals, nbins=nbins, closed=closed)
    hist_scene, hist_αs = plot_histogram(hist, cmap)

    hist_scene[Axis][:names, :axisnames] = (hname, "count")

    sc = AbstractPlotting.vbox(data_scene, hist_scene)

    selected_plot = setup_click(data_scene, 1)
    hist_idx = setup_click(hist_scene, 2)

    select_series(data_scene, selected_plot, data_αs, hist_αs, vals, hist, α, hα)
    select_bin(hist_idx, hist, hist_αs, data_αs, vals, closed=closed, α = α, hα = hα)
    display(sc)
    return data_scene, hist_scene
end

"""
    plot_histogram(hist, cmap) -> hist_scene, hist_αs
Plot a histogram where the transparency (α) of each bin can be changed
and return the scene together with the αs. The bins are colored
according to a colormap.
"""
function plot_histogram(hist::StatsBase.Histogram, cmap)
    c = to_colormap(cmap, length(hist.weights))
    hist_αs = [Observable(1.) for i in c]
    bincolor(αs...) = RGBAf0.(color.(c), αs)
    colors = lift(bincolor, hist_αs...)
    hist_scene = plot(hist, color=colors)
    return hist_scene, hist_αs
end

"""
    change_α(series_alpha, idxs, α = DEFAULT_α)
Given a vector of `Observable`s that represent the αs for some series, change
the elements with indices given by `idxs` to the value `α`.
This can be used to hide some series by using a low α value (default).
To restore the initial color, use `α = 1`.
"""
function change_α(series_alpha, idxs, α = DEFAULT_α)
    foreach(i -> series_alpha[i][] = α, idxs)
end

"""
    get_series_idx(selected_plot, scene)
Get the index of the `selected_plot` in `scene`.
"""
function get_series_idx(selected_plot, scene)
    # TODO: There is probably a better or more efficient way of doing this.
    plot_idx = findfirst(map(p->selected_plot === p, scene.plots))
    plot_idx
end

"""
    setup_click(scene, idx=1)
Given a `scene` return a `Observable` that listens to left clicks inside the scene.
The `idx` argument is used to index the tuple `(plt, click_idx)` which gives
the selected plot and the index of the selected element in the plot.
"""
function setup_click(scene, idx=1)
    selection = Observable{Any}(0)
    on(scene.events.mousebuttons) do buttons
        if ispressed(scene, Mouse.left) && AbstractPlotting.is_mouseinside(scene)
            plt, click_idx = AbstractPlotting.mouse_selection(scene)
            selection[] = (plt, click_idx)[idx]
        end
    end
    return selection
end

"""
    bin_with_val(val, hist)
Get the index of the bin in the histogram(`hist`) that contains the given value (`val`).
"""
bin_with_val(val, hist) = searchsortedfirst(hist.edges[1], val) - 1

"""
    idxs_in_bin(i, hist, val; closed=:left)
Given the values (`val`) which are histogramed in `hist`, find all the indices
which correspond to the values in the `i`-th bin.
"""
function idxs_in_bin(i, hist, val; closed=:left)
    h = hist.edges[1]
    inbin(x) = (closed == :left) ? h[i] ≤ x < h[i+1] : h[i] < x ≤ h[i+1]
    idx = findall(inbin, val)
    return idx
end

"""
    select_series(scene, selected_plot, data_αs, hist_αs, data, hist)
Setup selection of a series in a scatter plot and the corresponding histogram.
When a point of the scatter plot is clicked, the corresponding series is
highlighted (or selected) by changing the transparency of all the other series
(and corresponding histogram bins) to a very low value.
When the click is outside, the series is deselected, that is all the αs are
set back to 1.
"""
function select_series(scene, selected_plot, data_αs,
                       hist_αs, data, hist,
                       α = DEFAULT_α, hα = 0.2
                       )
    series_idx = map(get_series_idx, selected_plot, scene)
    on(series_idx) do i
        if !isa(i, Nothing)
            data_αs[i - 1][] = 1.0
            change_α(data_αs, setdiff(axes(data_αs, 1), i - 1), α)
            selected_bin = bin_with_val(data[i-1], hist)
            hist_αs[selected_bin][] = 1.0
            change_α(hist_αs, setdiff(axes(hist_αs, 1), selected_bin), hα)
        else
            change_α(data_αs, axes(data_αs, 1), 1.0)
            change_α(hist_αs, axes(hist_αs, 1), 1.0)
        end
        return nothing
    end
end

"""
    select_bin(hist_idx, hist, hist_αs, data_αs, data; closed=:left, α = DEFAULT_α)
Setup a selection of a histogram bin and the corresponding series in the
scatter plot. See also [`select_series`](@ref).
"""
function select_bin(hist_idx, hist, hist_αs, data_αs, data;
    closed=:left, α = DEFAULT_α, hα = 0.2)

    on(hist_idx) do i
        if i ≠ 0
            hist_αs[i][] = 1.0
            change_α(hist_αs, setdiff(axes(hist.weights, 1), i), hα)
            change_α(data_αs, idxs_in_bin(i, hist, data, closed=closed), 1.0)
            change_α(data_αs, setdiff(
                axes(data_αs, 1), idxs_in_bin(i, hist, data, closed=closed)
            ), α)
        else
            change_α(data_αs, axes(data_αs, 1), 1.0)
            change_α(hist_αs, axes(hist_αs, 1), 1.0)
        end
        return nothing
    end
end


text:
using Ramnet.Encoders

using StaticArrays

using ..Approximators: BinaryActionPolicy
using ..Buffers: MultiStepDynamicBuffer, add!

mutable struct FunctionalPG{OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}} <: AbstractAgent{O,Int}
  γ::Float64
  policy::BinaryActionPolicy{OS,T,O,E}
  buffer::MultiStepDynamicBuffer{O,Int}
  observation::Union{Nothing,O}
  done::Bool
  action::Union{Nothing,Int}
  min_obs::Union{Nothing,O}
  max_obs::Union{Nothing,O}
  rng::MersenneTwister
end

function FunctionalPG(::Type{O}, n, start_learning_rate, end_learning_rate, learning_rate_decay, epochs, discount, encoder::E; seed::Union{Nothing,UInt} = nothing) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
  !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
  rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

  n < 1 && throw(DomainError(n, "Tuple size must be at least 1"))

  policy = BinaryActionPolicy(O, n, encoder; start_learning_rate, end_learning_rate, learning_rate_decay, epochs, partitioner = :uniform_random, seed)

  buffer = MultiStepDynamicBuffer{O,Int}()

  FunctionalPG{OS,T,O,E}(discount, policy, buffer, nothing, false, nothing, nothing, nothing, rng)
end

function FunctionalPG(
  env, encoder::E; tuple_size, start_learning_rate, end_learning_rate, learning_rate_decay, epochs = 1, discount = 1, seed = nothing) where {T<:Real,E<:AbstractEncoder{T}}

  FunctionalPG(observation_type(env), tuple_size, start_learning_rate, end_learning_rate, learning_rate_decay, epochs, discount, encoder; seed)
end

function observe!(agent::FunctionalPG{OS,T,O,E}, observation::O) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
  agent.observation = observation
  agent.done = false

  # Finding the edges of the domain
  if isnothing(agent.min_obs)
    agent.min_obs = observation
  else
    agent.min_obs = min.(agent.min_obs, observation)
  end

  if isnothing(agent.max_obs)
    agent.max_obs = observation
  else
    agent.max_obs = max.(agent.max_obs, observation)
  end

  nothing
end

# TODO: All observe! methods are the same for all agents. Generalize.
# TODO: This method does not need to take in the action
function observe!(agent::FunctionalPG{OS,T,O,E}, action::Int, reward::Float64, observation::O, done::Bool) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
  add!(agent.buffer, agent.observation, agent.action, reward)

  agent.observation = observation
  agent.done = done

  # Finding the edges of the domain
  agent.min_obs = min.(agent.min_obs)
  agent.max_obs = max.(agent.max_obs)

  nothing
end

function select_action(agent::FunctionalPG{OS,T,O,E}, observation::O; deterministic::Bool = false) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
  return select_action(agent.policy, observation; deterministic)
end

function select_action!(agent::FunctionalPG{OS,T,O,E}, observation::O) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
  agent.action = select_action(agent, observation)

  return agent.action
end

function update!(agent::FunctionalPG{OS,T,O,E}) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
  if agent.done
    for _ in 1:agent.policy.opt.epochs
      G = 0.0
      for transition in Iterators.reverse(agent.buffer)
        G = transition.reward + agent.γ * G

        update!(agent.policy, transition.observation, transition.action, G)
      end
    end

    reset!(agent.buffer)
  end

  nothing
end

function Agents.reset!(agent::FunctionalPG{OS,T,O,E}; seed::Union{Nothing,UInt} = nothing) where {OS,T<:Real,O<:StaticArray{Tuple{OS},T,1},E<:AbstractEncoder{T}}
  reset!(agent.policy; seed)
  reset!(agent.buffer)

  agent.observation = nothing
  agent.action = nothing
  agent.done = false

  nothing
end


text:
#tbeg = time()
using LinearAlgebra
BLAS.set_num_threads(4)
using Flux, Statistics # Flux.Data.MNIST
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated, partition
using BSON, HDF5, JLD, Random
using MLDataUtils
using NPZ
include("genDataScripts.jl")
n_sites = 100
periods = [2,3,4,5,6,7,8,9,10,11]
all_signals = gen_signals(n_sites, periods)

neu = 100
#trainingSizes = [50, 100, 250, 500, 1000, 2000, 4000]
trainSizes = [50, 100, 250, 500, 1000, 2000, 4000]
sigmasAll = collect(0.05:0.05:4)
indice = parse(Int64, ARGS[1])
#σv = round(σvt*0.05, digits=3)
#σv = 2.5 # 1.0

grid = []
for ts in trainSizes, σv in sigmasAll
    push!(grid, [ts, σv])
end
tssize = Int64(grid[indice][1])
trainingSizes = [tssize]
σv = grid[indice][2]
iteraciones = 50

# HHG DATA
x_data9k = npzread("./data/data_simulation2_optimized_9669.npy")
labels9k = npzread("./data/labels_simulation2_optimized_9669.npy")
ycold9k = Int64.(labels9k[:,1,1])
nsamples = size(x_data9k)[1]
npts = 100
data_log9k = zeros(npts, nsamples)
for i in 1:nsamples
    input_log = log10.(x_data9k[i,:,2])
    input_log = input_log .+ minimum(input_log)*(-1)
    #input_log = input_log/maximum(input_log)
    data_log9k[:,i] .= input_log
end
accuracy_tcold(x,y,model) = mean(onecold(model(x)) .== y)
function signs100k(σv)
    Random.seed!(1234)
    signals_mas_ruido = zeros(npts,nsamples)
    for i in 1:nsamples
        signals_mas_ruido[:,i] .= signal_plus_noise(data_log9k[:,i], σv, true)
    end
    for _ in 1:9
        signals_mas_ruidox = zeros(npts, nsamples)
        for i in 1:nsamples
            signals_mas_ruidox[:,i] .= signal_plus_noise(data_log9k[:,i], σv, true)
        end
        signals_mas_ruido = hcat(signals_mas_ruido, signals_mas_ruidox)
    end
    signals_mas_ruido
end
ycold9kx = repeat(ycold9k, 10)

# Synthetic data
function synthData(σv)
    all_signals = gen_signals(n_sites, periods)
    Random.seed!(1234)
    num_test = 10000
    x_test, y_test = gen_data_set(all_signals, σv, num_test,true)
    xtest_temp, ytest = reshape(x_test, (100, num_test, 10)), reshape(y_test, (num_test, 10))
    x_test = reshape(xtest_temp, (n_sites, 10*num_test))
    labels_test = reshape(ytest, (10*num_test))
    x_test, labels_test
end

Random.seed!(1234)
t_size = 1_000_000 # 4_000_000 old
num_train = Int64(round.(t_size *0.8/10)) # number of samples for each class,
#this just works for 10 classes
num_val = Int64(round.(t_size *0.2/10))
num_test = 1000

x_test, y_test = gen_data_set(all_signals, σv, num_test,true)
x_train, y_train = gen_data_set(all_signals, σv, num_train,true)
x_vali, y_vali = gen_data_set(all_signals, σv, num_val,true)

x_test = reshape(x_test, (100, num_test, 10))
y_test = reshape(y_test, (num_test, 10))
x_train = reshape(x_train, (100, num_train, 10))
y_train = reshape(y_train, (num_train, 10))
x_vali = reshape(x_vali, (100, num_val, 10))
y_vali = reshape(y_vali, (num_val, 10))

xtest_temp = x_test[:,1:num_test,:]
ytest_temp = y_test[1:num_test,:]
x_test = reshape(xtest_temp, (n_sites, 10*num_test))
labels_test = reshape(ytest_temp, (10*num_test))

X_test = x_test # |> gpu
Y_test = onehotbatch(labels_test, 1:10); #|> gpu
for t_size in trainingSizes
    tmtrain = time()
    num_train = Int64(round.(t_size *0.8/10)) #number of samples for each class,
    #this just works for 10 classes
    num_val = Int64(round.(t_size *0.2/10))

    xtrain_temp = x_train[:,1:num_train,:]
    ytrain_temp = y_train[1:num_train,:]

    x_data = reshape(xtrain_temp, (n_sites, 10*num_train))
    labels = reshape(ytrain_temp, (10*num_train))

    xvali_temp = x_vali[:,1:num_val,:]
    yvali_temp = y_vali[1:num_val,:]

    x_val = reshape(xvali_temp, (n_sites, 10*num_val))
    labels_val = reshape(yvali_temp, (10*num_val))

    # Data preparation
    X = x_data # |> gpu
    Y = onehotbatch(labels, 1:10); #|> gpu
    #dataset = [(X, Y)]
    X_val = x_val # |> gpu
    Y_val = onehotbatch(labels_val, 1:10); #|> gpu

    #mini-batch
    batch_size = 256 #512 is ok too, similar performance at test accuracy
    mb_idxs = partition(1:size(X)[2], batch_size)
    #train_set = [(X[:,p], Y[:,p]) for p in mb_idxs]
    for itera in 51:iteraciones+50
        # saving training history
        epochs = 1000
        t_loss = fill(NaN, epochs)
        v_loss = fill(NaN, epochs)
        va_acc = fill(NaN, epochs)

        #saving best models
        dict_models = Dict()

        #dict_models["model_iter"] = NaN
        dict_models["iter_epoch_num"] = 0

        #saving best model
        best_val = 5.0
        last_improvement = 0
        patience=50

        Random.seed!(10*itera)

        model = Chain(
          Dense(n_sites, neu, Flux.relu),
          Dense(neu, 10),
          Flux.softmax) # |> gpu
        opt = ADAM()
        loss(x, y) = crossentropy(model(x), y)
        global accuracyT(x, y) = mean(onecold(model(x)) .== onecold(y))

        for epoch_indx in 1:epochs
            xs, ys = shuffleobs((X, Y))
            train_set = [(xs[:,p], ys[:,p]) for p in mb_idxs]
            Flux.train!(loss, params(model), train_set, opt)

            validation_loss = Tracker.data(loss(X_val, Y_val))
            t_loss[epoch_indx] = Tracker.data(loss(X, Y))
            v_loss[epoch_indx] = validation_loss
            va_acc[epoch_indx] = accuracyT(X_val, Y_val)

            # If this is the best val_loss we've seen so far, save the model out
            if validation_loss <= best_val
                #@info(" -> New best val_loss! Saving model out to model_iter$(r).bson")
                last_improvement = epoch_indx
                dict_models["iter_epoch_num"] = last_improvement
                best_val = validation_loss
            end
            if epoch_indx - last_improvement >= patience
                #@info(" -> Early-exiting iteration $(r) and epoch $(epoch_indx): no more patience")
                break
            end
        end
        test_acc = accuracyT(X_test, Y_test)

        BSON.@save "modelsN$(neu)ESP_ts$(t_size)_σ$(σv)_ite$(itera).bson" model t_loss v_loss va_acc test_acc

        error_TS_ruido = zeros(length(sigmasAll))
        for (i,σvd) in enumerate(sigmasAll)
            signals_mas_ruido = signs100k(σvd)
            error_TS_ruido[i] = (1.0 .- accuracy_tcold(signals_mas_ruido, ycold9kx, model))
        end
        npzwrite("error_hhg_$(t_size)_1n_sig_$(σv)_ite$(itera).npy", error_TS_ruido)

        error_TS_ruidox = zeros(length(sigmasAll))
        for (i,σvd) in enumerate(sigmasAll)
            x_test, labels_test = synthData(σvd)
            error_TS_ruidox[i] = (1.0 .- accuracy_tcold(x_test, labels_test, model))
        end
        npzwrite("error_GtoG100k_$(t_size)_1n_sig_$(σv)_ite$(itera).npy", error_TS_ruidox)
    end
end
#println(time() - tbeg)

max_stars_repo_name: lazarusA/noisySignals

max_issues_repo_name: lazarusA/noisySignals

max_forks_repo_name: lazarusA/noisySignals


text:
"""
    gradient_hessian_equal_rates(λ, i, j, t)

Compute the gradient of the log-probability of a simple birth and death process
evaluated at the point ``η = (λ, λ)^{\\prime}``. Variable `i` is
the initial population size, `j` is the final population size, and `t` is the
elapsed time.
"""
function gradient_hessian_equal_rates(
  λ::F,
  i::F,
  j::F,
  t::F
)::Tuple{Vector{F}, Symmetric{F, Matrix{F}}} where {
  F <: AbstractFloat
}
  r1 = zero(F)
  r2 = zero(F)
  r3 = zero(F)
  r4 = zero(F)
  r5 = zero(F)

  θ = λ * t
  y = θ / (1 + θ)

  if j > 1
    if i > 1
      a, b = (j <= i) ? (i, j) : (j, i)

      q1, q2, q3 = if θ <= 1
        log_hypergeometric_joint(a, b, logexpm1(-2 * log(θ)))
      else
        log_meixner_ortho_poly_joint(a, b, 2 * log(θ))
      end

      w1 = i * j * exp(q2 - q1) / (i + j - 1)
      w2 = (i - 1) * (j - 1) * exp(q3 - q2) / (i + j - 2)
      w3 = λ * θ^2

      x1 = (i + j) * y / 2

      # firt-order partial derivative with respect to λ
      r1 = -(w1 + θ^2 * (x1 - j)) / w3

      # first-order partial derivative with respect to μ
      r2 = -(w1 + θ^2 * (x1 - i)) / w3

      x1 = (i + j) * y^2 / 12
      x2 = w1 * (w1 - w2 + θ^2 * (θ^2 - 12) / 6)

      c1 = 2 * θ - 1
      c2 = 2 * θ + 5

      # second-order partial derivative with respect to λ twice
      r3 = -(x2 + θ^4 * (x1 * c1 + j)) / w3^2

      # second-order partial derivative with respect to λ and then μ
      r4 = -(w1 * (w1 - w2 - θ^2 * (θ^2 + 6) / 6) - θ^4 * x1 * c2) / w3^2

      # second-order partial derivative with respect to μ twice
      r5 = -(x2 + θ^4 * (x1 * c1 + i)) / w3^2
    else
      r1 = (j - 1) / (λ * (1 + θ)) + ((j - 3) / (2 * λ)) * y
      r2 = -((j + 1) / (2 * λ)) * y
      r3 = -(12 * (j - 1) +
             θ * (24 * (j - 1) +
                  θ * (11 * (j - 1) +
                       θ * (2 * (j + 3) +
                            θ * 2)))) / (12 * (λ * (1 + θ))^2)
      r4 = (5 * j + 7 + 2 * θ * (j + 3 + θ)) * y^2 / (12 * λ^2)
      r5 = (j - 1 - 2 * θ * (j + 3 + θ)) * y^2 / (12 * λ^2)
    end
  elseif j == 0
    a1 = i / (2 * λ)
    a2 = i / (12 * λ^2)

    r1 = -a1 * y
    r2 = a1 * ((2 + θ) / (1 + θ))
    r3 = a2 * (1 - 2 * θ) * y^2
    r4 = a2 * (5 + 2 * θ) * y^2
    r5 = -a2 * (12 + θ * (24 + θ * (11 + θ * 2))) / (1 + θ)^2
  elseif j == 1
    r1 = -((i + 1) / (2 * λ)) * y
    r2 = (i - 1) / (λ * (1 + θ)) + ((i - 3) / (2 * λ)) * y
    r3 = (i - 1 - 2 * θ * (i + 3 + θ)) * y^2 / (12 * λ^2)
    r4 = (5 * i + 7 + 2 * θ * (i + 3 + θ)) * y^2 / (12 * λ^2)
    r5 = - (12 * (i - 1) +
            θ * (24 * (i - 1) +
                 θ * (11 * (i - 1) +
                      θ * (2 * (i + 3) +
                           θ * 2)))) / (12 * (λ * (1 + θ))^2)
  end

  [r1, r2], Symmetric([r3 r4; r4 r5])
end


text:
module Effects

using Base.Iterators: product
using DataFrames
using LinearAlgebra
using Statistics
using StatsModels
using StatsBase
using Tables

using StatsModels: AbstractTerm

export effects
export effects!

include("typical.jl")
include("regressionmodel.jl")

end # module


text:
using Bitly
using Test

const testurl = "https://docs.julialang.org/en/v1/"

@testset "Shorten" begin
    b = BitlyToken()
    @test b.url == Bitly.BITLY_API_BASE
    S = shorten(b,testurl)

    @test S.link == "https://bit.ly/2Us5Vl7"

    @test S.response["long_url"] == testurl
end

@testset "Expand" begin
    b = BitlyToken()
    S = expand(b,"https://bit.ly/2Us5Vl7")
    @test S.long_url == testurl
end


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

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

# ╔═╡ 21b1f4b6-2172-4292-8249-35a615db4f59
begin
    using Pkg
	Pkg.activate(mktempdir())
	# Pkg.add("Revise");using Revise
	Pkg.add("PlutoUI")
	using PlutoUI
	Pkg.add(name="PlutoCanvasPlot",version="0.0.4")
#	Pkg.develop("PlutoCanvasPlot")
	using PlutoCanvasPlot

end

# ╔═╡ c54da25c-0dbf-41d1-bd7f-657b115bf6bf
md"""
Test Notebook for [PlutoCanvasPlot](https://github.com/j-fu/PlutoCanvasPlot.jl)

So far, this package is in an early state. The current version  of the code is available via the registry `https://github.com/j-fu/PackageNursery.git` .

So, in order to run this notebook, you need to add this registry to you Julia environment.

```
pkg> registry add https://github.com/j-fu/PackageNursery.jl.git
```

This step can be safely undone by removing `.julia/registries/PackageNursery` in your
Julia folder.

"""

# ╔═╡ 825f6068-0f02-44e4-b083-1be3eb4c764d
let
	X=0:0.1:10
	p=CanvasPlot(;xrange=X,yrange=-1:1,resolution=(600,200))
	polyline!(p,X,sin.(X))
	polyline!(p,X,cos.(X))
    p
end

# ╔═╡ 2a5510d3-7355-4095-ab43-ad6d75bba090
let
	p=CanvasPlot(resolution=(600,300))

	axis!(p; xtics=(0:2.5:10),ytics=(-1:0.5:1))
	X=0:0.1:10
	linecolor!(p,1,0,0)
   	polyline!(p,X,sin.(X))
	linecolor!(p,0,1,0)
	polyline!(p,X,cos.(X))
end

# ╔═╡ 405b81de-3fd9-48f7-81a7-dc8f06dfda9c
let
	p=CanvasPlot(resolution=(300,300),xrange=-1:1,yrange=-1:1)
	T=0:0.01:2π
	X=0.45.*sin.(T)
	Y=0.45.*cos.(T)
	fillcolor!(p,0,1,0)
	polygon!(p,X,Y)
end

# ╔═╡ 83c16003-8196-4f2c-a602-c298f27475fb
begin
	@bind clock_t Clock(interval=0.01)
end

# ╔═╡ fd7c959b-94d6-43ff-b535-49b0883874d5
let
	p=CanvasPlot(resolution=(300,300),xrange=-1:1,yrange=-1:1)
	time=0.005*clock_t
	dt=0.01
    T=time:dt:2π/2+time
	px=5+sin(time/10)
	py=6+cos(time/10)
	X=sin.(px .*T)
	Y=cos.(py .*T)
	polyline!(p,X,Y)
end

# ╔═╡ Cell order:
# ╟─c54da25c-0dbf-41d1-bd7f-657b115bf6bf
# ╠═21b1f4b6-2172-4292-8249-35a615db4f59
# ╠═825f6068-0f02-44e4-b083-1be3eb4c764d
# ╠═2a5510d3-7355-4095-ab43-ad6d75bba090
# ╠═405b81de-3fd9-48f7-81a7-dc8f06dfda9c
# ╟─83c16003-8196-4f2c-a602-c298f27475fb
# ╠═fd7c959b-94d6-43ff-b535-49b0883874d5


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using SymPy

# %%
function minpoly(x, y, f, g, h, N)
    m = sympy.Poly(f, x).degree()
    n = sympy.Poly(g, y).degree()
    K = sympy.Poly(h, x).degree()
    L = sympy.Poly(h, y).degree()
    λ = symbols("λ[0:$N]")
    z = symbols("z")
    P = sum(λ[begin+k]*z^k for k in 0:N-1) + z^N
    p = P(z => h).expand()
    @time for k in K*N:-1:m
        p = p(x^k => x^k - x^(k-m)*f).expand()
    end
    @time for l in L*N:-1:n
        p = p(y^l => y^l - y^(l-n)*g).expand()
    end
    Q = sympy.Poly(p, x, y)
    C = Q.coeffs()
    @time sol = solve(C, λ)
    P = sum(get(sol, λ[begin+k], λ[begin+k])*z^k for k in 0:N-1) + z^N
    [
        sympy.Poly(f, x)
        sympy.Poly(g, y)
        sympy.Poly(h, x, y)
        sympy.Poly(P, z)
    ]
end

# %%
@vars a p x y
f = x^2 - a
g = y^2 - p
h = x + y
minpoly(x, y, f, g, h, 4)

# %%
@vars a p x y
f = x^2 - a
g = y^2 - p
h = x * y
minpoly(x, y, f, g, h, 2)

# %%
@vars a p x y
f = x^2 - a
g = y^3 - p
h = x + y
minpoly(x, y, f, g, h, 6)

# %%
@vars a p x y
f = x^2 - a
g = y^3 - p
h = x * y
minpoly(x, y, f, g, h, 6)

# %%
@vars a p q x y
f = x^2 - a
g = y^3 + p*y + q
h = x + y
minpoly(x, y, f, g, h, 6)

# %%
@vars a p q x y
f = x^2 - a
g = y^3 + p*y + q
h = x * y
minpoly(x, y, f, g, h, 6)

# %%
@vars a p q r x y
f = x^2 - a
g = y^4 + p*y^2 + q*y + r
h = x + y
minpoly(x, y, f, g, h, 8)

# %%
@vars a p q r x y
f = x^2 - a
g = y^4 + p*y^2 + q*y + r
h = x * y
minpoly(x, y, f, g, h, 8)

# %%
@vars a p q r s x y
f = x^2 - a
g = y^5 + p*y^3 + q*y^2 + r*y + s
h = x + y
minpoly(x, y, f, g, h, 10)

# %%
@vars a p q r s x y
f = x^2 - a
g = y^5 + p*y^3 + q*y^2 + r*y + s
h = (x * y)^2
minpoly(x, y, f, g, h, 5)

# %%
@vars a p x y
f = x^3 - a
g = y^3 - p
h = (x + y)^3
minpoly(x, y, f, g, h, 3)

# %%
@vars a p x y
f = x^3 - a
g = y^3 - p
h = x * y
minpoly(x, y, f, g, h, 3)

# %%
@vars a b p q x y
f = x^3 + a*x + b
g = y^3 + p*y + q
h = x + y
minpoly(x, y, f, g, h, 9)

# %%
@vars a b p q x y
f = x^3 + a*x + b
g = y^3 + p*y + q
h = x * y
minpoly(x, y, f, g, h, 9)

# %%
@vars a p q x y
f = x^3 - a
g = y^4 + p*y + q
h = x + y
minpoly(x, y, f, g, h, 12)

# %%
@vars a p q x y
f = x^3 - a
g = y^4 + p*y + q
h = x * y
minpoly(x, y, f, g, h, 12)

# %%

max_stars_repo_path: 0023/minimal polynomials of sums and products.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
"""
    fixedpoint_gen(func, u0, p, ranges; MaxIter=1000, UniqSpac=8)
    fixedpoint_gen(func, p, ranges; MaxIter=1000, UniqSpac=8)

Find steady states and categorized stability with Jacobian.

Arguements
----------
- `func`{Function}: An ODE function
- `u0`{Vector}: Initial Values (used for creating `DEsteady`)
- `p`{Vector}: Parameters
- `ranges`{Vector of range}: define the domain and grid for Searching
- `MaxIter`: Maximum iteration for solving the root
- `UniqSpac`{Int}: Rounding the solution.
"""
function fixedpoint_gen(func, u0, p, ranges; MaxIter=10000, UniqSpac=8)

ssmethod = SSRootfind(nlsolve = (f,u0,abstol) -> (
    res=nlsolve(f,u0, autodiff=:forward, iterations= MaxIter,ftol=abstol);
    res.zero
    ))

# Define a problem
de = DEsteady(func=func, p=p, u0= u0, method= ssmethod)

j_gen = jacobian(de) # jacobian generator

# Searching method and domain
param_gen = ParameterGrid(ranges)

# Solve
sols = solve(de, param_gen)

# Remove similar solutions
steadies = unique(sols; tol_digit=UniqSpac)

# Jacobian
jac_ms = j_gen.(steadies)

# Stability
stab_modes = StabilityType.(jac_ms)

return steadies, jac_ms, stab_modes

end

function fixedpoint_gen(func, p, ranges; kwag...)
    u0 = [i[1] for i in ranges]
    return fixedpoint_gen(func, u0, p, ranges; kwag...)
end


text:
include("src/base.jl")
include("src/envs/toytext_mdp.jl")
include("src/policies/value_table_policy.jl")
include("src/solvers/ev_sarsa_policy_solver.jl")

gym_env = ToyTextMDP(gym.make("Taxi-v2"))
solver = EVSarsaPolicySolver(5000; max_frame_iterations = 100, epsilon = 0.25, epsilon_discount = 0.995, learning_rate = 0.1, discount = 0.99)
policy = solve(solver, gym_env, verbose = true)

# println(policy.value_map)

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
"""
    set_properties!(dict, obj; opts...)
    set!(...)

Set properties of an object `obj` given options (`opts`) of the form
`optname=value`, and applying it through an appropriate application function
stored in the dictionary `dict`.
Note that the dictionary of property-setting-functions also contains
"pre-conditioners" which are functions that check that the values passed to
properties are sensible and convert them to sensible types if relevant.
This reduces code duplication and allows to reduce specialization.

Note: `set_properties!` (to avoid ambiguities in the code) is used internally,
    `set!` is exported (less verbose and not ambiguous).
"""
function set_properties!(
            dict::LittleDict{Symbol, Pair{Function, Function}},
            obj;
            opts...
        )::Nothing

    for optname in keys(opts)
        optvalue = opts[optname]
        argcheck, setprop! = get(dict, get(PROP_ALIAS, optname, optname)) do
            throw(
                UnknownOptionError(optname, obj)
            )
        end
        setprop!(obj, argcheck(optvalue, optname))
    end
    return
end
set!(obj; opts...) = set_properties!(obj; opts...)

####
#### Value checkers for set_properties functions the symbol corresponds to the name
#### of the option that is being modified
####

# No Check
I(x, ::Symbol) = x

# Bool
bool(x::Bool,      ::Symbol) = x
bool(x::Int,       ::Symbol) = !iszero(x)
bool(x::Listable, o::Symbol) = b.(x, o)
bool(x,           o::Symbol) = throw(
    OptionValueError(o, x, "((List of) Bool)")
)

# Reversal of input (e.g. legend box/nobox)
not(x::Bool,  ::Symbol) = !x
not(x,       o::Symbol) = throw(
    OptionValueError(o, x, "((List of) Bool)")
)

# string or vector
str(x::String,    ::Symbol)    = x
str(x::Listable, o::Symbol) = str.(x, o)
str(x,           o::Symbol) = throw(
    OptionValueError(o, x, "((List of) String)")
)

# Lowercase string or strings
lc(x::String,    ::Symbol) = lowercase(x)
lc(x::Listable, o::Symbol) = lc.(x, o)
lc(x,           o::Symbol) = throw(
    OptionValueError(o, x, "((List of) String)")
)

# Cast to float
fl(x::Real,      ::Symbol) = Float64(x)
fl(x::Listable, o::Symbol) = fl.(x, o)
fl(x,           o::Symbol) = throw(
    OptionValueError(o, x, "((List of) Real)")
)

pos(x::Real, o::Symbol) = (x < 0 && throw(
    OptionValueError(o, x, "(positive)")
);)

# Check positive and cast to float
posfl(x::Real,     o::Symbol) = (pos(x, o); Float64(x))
posfl(x::Listable, o::Symbol) = posfl.(x, o)
posfl(x,           o::Symbol) = throw(
    OptionValueError(o, x, "((List of) positive Real)")
)

# check positive + integer
posint(x::Int, o::Symbol) = (pos(x, o); x)
posint(x,      o::Symbol) = throw(
    OptionValueError(o, x, "(positive integer)")
)

# color
col(c::Color,     ::Symbol) = c
col(s::String,    ::Symbol) = parse(Color, s)
col(v::Listable, o::Symbol) = col.(v, o)
col(x,           o::Symbol) = throw(
    OptionValueError(o, x, "((List of) Color or String-Color)")
)

"""
    opcol(x, o)

Internal function to process an optional color, i.e. `"none"` can be passed
resulting in `nothing` being passed.
This is relevant for instance when setting the background color of a figure
to "none" (transparent figure). If `"none"` is passed, the current figure
is checked for transparency properties if set to false, a warning will be
raised and the color will be defaulted to `"white"` if the transparency
setting is unset, it will be set to `true`.
"""
function opcol(x::String, o::Symbol)::Option{Color}
    if lowercase(x) == "none"
        isdef(gcf().transparency) || (gcf().transparency = true)
        if !gcf().transparency
            @warn """
                Transparent background is only supported when the figure
                has its transparency property set to 'true'.
                """
            return col("white", o) # fully opaque
        end
        return nothing
    end
    return col(x, o)
end
opcol(c, o) = col(c, o)

"""
    alpha(α, s)

Internal function to process an alpha parameter (transparency). The current
figure is checked for transparency properties, if set to false, a warning will
be raised and the `α` will be ignored (fully opaque). If the transparency
setting is unset, it will be set to `true`.
Accepted values are strictly between `0` and `1`.
For completely transparent, see using `"none"` in the color description
(for instance `Figure(bgcol="none")`).
"""
function alpha(α::Float64, optname::Symbol)::Float64
    isdef(gcf().transparency) || (gcf().transparency=true)
    if !gcf().transparency
        @warn """
            Transparent colors are only supported when the figure
            has its transparency property set to 'true'. Ignoring α.
            """
        return 1.0 # fully opaque
    end
    0 < α <= 1 || throw(
        OptionValueError("alpha", α, "(number in (0, 1])")
    )
    return α
end
alpha(α::Real, o::Symbol) = alpha(fl(α), o)
alpha(α,       o::Symbol) = throw(
    OptionValueError("alpha", α, "(number in (0, 1])")
)

####
#### Pickers for GROUPED objects
#### --> gline, gbar, scatter, boxplot, bar
####

# this behaves like ∘ except it splats the output of g, this is useful in
# set_properties when wanting to use a set function but it has to be applied
# on a subfield of the object, then we apply a "selector" first to pick that
# subfield before applying the setter.
# The effect would be something like setter(picker(obj), value) and is
# written (setter ⊙ picker)
⊙(f::Function, g::Function) = (args...) -> f(g(args...)...)

###############################################################
####
#### Options for STYLE
####
###############################################################

pick_lstyles(o, v) = (o.linestyles,   v)
pick_mstyles(o, v) = (o.markerstyles, v)
pick_bstyles(o, v) = (o.barstyles,    v)
pick_blstyle(b, v) = (b.boxstyles,    v, :blstyle)
pick_mlstyle(b, v) = (b.boxstyles,    v, :mlstyle)
pick_mmstyle(b, v) = (b.boxstyles,    v, :mmstyle)
pick_omstyle(b, v) = (b.boxstyles,    v, :omstyle)

# style.jl
# --------
# set_color! <- color (option for Figure/Legend)
# set_font! <- string
# set_hei!  <- F64, positive

const PROP_ALIAS = LittleDict(
    :color              => :col,
    :colour             => :col,
    :colors             => :col,
    :colours            => :col,
    #
    :textcol            => :tc,
    :textcolor          => :tc,
    :textcolour         => :tc,
    #
    :lc                 => :col,
    :linecol            => :col,
    :linecolor          => :col,
    :linecolors         => :col,
    :linecolour         => :col,
    :linecolours        => :col,
    #
    :mcol               => :mc,
    :markercol          => :mc,
    :markercolor        => :mc,
    :markercolour       => :mc,
    #
    :ecol               => :ec,
    :edgecol            => :ec,
    :edgecolor          => :ec,
    :edgecolour         => :ec,
    :edgecolors         => :ec,
    :edgecolours        => :ec,
    #
    :fill               => :col,
    :fills              => :col,
    #
    :lstyle             => :ls,
    :linestyle          => :ls,
    :lstyles            => :ls,
    :linestyles         => :ls,
    #
    :lwidth             => :lw,
    :linewidth          => :lw,
    :lwidths            => :lw,
    :linewidths         => :lw,
    #
    :smooths            => :smooth,
    #
    :markers            => :marker,
    #
    :msize              => :ms,
    :markersize         => :ms,
    :msizes             => :ms,
    :markersizes        => :ms,
    #
    :bwidth             => :width,
    :barwidth           => :width,
    #
    :position           => :pos,
    #
    :bgcolor            => :bgcol,
    :bgcolour           => :bgcol,
    :background         => :bgcol,
    #
    :length             => :len,
    #
    :symticks           => :sym,
    #
    :distance           => :dist,
    #
    :tickscol           => :tcol,
    :tickscolor         => :tcol,
    :tickscolour        => :tcol,
    #
    :key                => :label,
    :keys               => :label,
    :labels             => :label,
    #
    :from               => :xmin,
    :to                 => :xmax,
    #
    :norm               => :scaling,
    :nbins              => :bins,
    :horizontal         => :horiz,
    #
    :box_width          => :bw,
    :box_widths         => :bw,
    #
    :whisker_width      => :ww,
    :whisker_widths     => :ww,
    :box_wwidth         => :ww,
    :box_wwidths        => :ww,
    :box_whisker_width  => :ww,
    :box_whisker_widths => :ww,
    #
    :whisker            => :wl,
    :whisker_length     => :wl,
    :whiskers_length    => :wl,
    :whisker_lengths    => :wl,
    :whiskers_lengths   => :wl,
    #
    :box_ls             => :bls,
    :box_lstyle         => :bls,
    :box_lstyles        => :bls,
    :box_linestyle      => :bls,
    :box_linestyles     => :bls,
    #
    :box_lw             => :blw,
    :box_lwidth         => :blw,
    :box_lwidths        => :blw,
    :box_linewidths     => :blw,
    #
    :box_col            => :bc,
    :box_cols           => :bc,
    :box_color          => :bc,
    :box_colour         => :bc,
    :box_colors         => :bc,
    :box_colours        => :bc,
    #
    :med_ls             => :medls,
    :med_lstyle         => :medls,
    :med_lstyles        => :medls,
    :med_linestyle      => :medls,
    :med_linestyles     => :medls,
    #
    :med_lw             => :medlw,
    :med_lwidth         => :medlw,
    :med_lwidths        => :medlw,
    :med_linewidth      => :medlw,
    :med_linewidths     => :medlw,
    #
    :med_col            => :medc,
    :med_cols           => :medc,
    :med_color          => :medc,
    :med_colors         => :medc,
    #
    :mean_marker        => :mm,
    :mean_markers       => :mm,
    #
    :mean_msize         => :mms,
    :mean_markersize    => :mms,
    :mean_markersizes   => :mms,
    #
    :mean_mcol          => :mmc,
    :mean_mcols         => :mmc,
    :mean_markercol     => :mmc,
    :mean_markercolor   => :mmc,
    :mean_markercolors  => :mmc,
    :mean_markercolour  => :mmc,
    :mean_markercolours => :mmc,
    #
    :out_marker         => :om,
    :out_markers        => :om,
    #
    :out_msize          => :oms,
    :out_msizes         => :oms,
    :out_markersize     => :oms,
    :out_markersizes    => :oms,
    #
    :out_mcol           => :omc,
    :out_mcols          => :omc,
    :out_markercol      => :omc,
    :out_markercolor    => :omc,
    :out_markercolors   => :omc,
    :out_markercolour   => :omc,
    :out_markercolours  => :omc,
    #
    :colormap           => :cmap,
    :colourmap          => :cmap,
    #
    :res                => :pixels,
    :resolution         => :pixels,
    #
    :hastex => :tex,
    :latex => :tex,
    :haslatex => :tex,
    #
    :transparency => :alpha,
    :transparent  => :alpha,
    #
    :texpreamble => :preamble,
)

const OptsDict = LittleDict{Symbol,Pair{Function,Function}}

const TEXTSTYLE_OPTS = OptsDict(
    :font     => lc    => set_font!,
    :fontsize => posfl => set_hei!,
    :tc       => col   => set_textcolor!,
    )

const LINESTYLE_OPTS = OptsDict(
    :ls     => I     => set_lstyle!,
    :lw     => posfl => set_lwidth!,
    :smooth => bool  => set_smooth!,
    :col    => col   => set_color!,
    )

const GLINESTYLE_OPTS = OptsDict(
    :ls     => I     => set_lstyles! ⊙ pick_lstyles,
    :lw     => posfl => set_lwidths! ⊙ pick_lstyles,
    :smooth => bool  => set_smooths! ⊙ pick_lstyles,
    :col    => col   => set_colors!  ⊙ pick_lstyles,
    )

const MARKERSTYLE_OPTS = OptsDict(
    :marker => lc    => set_marker!,
    :ms     => posfl => set_msize!,
    :mc     => col   => set_mcol!,
    )

const GMARKERSTYLE_OPTS = OptsDict(
    :marker => lc    => set_markers! ⊙ pick_mstyles,
    :ms     => posfl => set_msizes!  ⊙ pick_mstyles,
    :mc     => col   => set_mcols!   ⊙ pick_mstyles,
    )

const BARSTYLE_OPTS = OptsDict(
    :ec => col => set_color!, # color
    )

const GBARSTYLE_OPTS = OptsDict(
    :col   => col   => set_fills!  ⊙ pick_bstyles,
    :ec    => col   => set_colors! ⊙ pick_bstyles,
    :width => posfl => set_bwidth!,
    )

const FILLSTYLE_OPTS = OptsDict(
    :col   => col   => set_fill!,
    :alpha => alpha => set_alpha!,
    )

###############################################################
####
#### Options for AX_ELEMS
####
###############################################################

const TITLE_OPTS = OptsDict(
    :dist => posfl => set_dist!
    )
merge!(TITLE_OPTS, TEXTSTYLE_OPTS)
set_properties!(t::Title; opts...) = set_properties!(TITLE_OPTS, t; opts...)

const LEGEND_OPTS = OptsDict(
    # position of the legend
    :pos     => lc    => set_position!, # set_legend
    # toggle legend on/off
    :off     => bool  => set_off!,      # .
    :on      => not   => set_off!,      # .
    # toggle surrounding box on/off
    :nobox   => bool  => set_nobox!,    # .
    :box     => not   => set_nobox!,    # .
    # specify margins (internal) and offset (external)
    :margins => fl    => set_margins!,  # .
    :offset  => fl    => set_offset!,   # .
    # background color and alpha
    :bgcol   => opcol => set_color!,    # set_style
    :bgalpha => alpha => set_alpha!,    # .
    )
merge!(LEGEND_OPTS, TEXTSTYLE_OPTS)
set_properties!(l::Legend; opts...) = set_properties!(LEGEND_OPTS, l; opts...)

const TICKS_OPTS = OptsDict(
    # ticks
    :off        => bool => set_off!,
    :on         => not  => set_off!,
    :len        => fl   => set_length!,
    :sym        => bool => set_symticks!,
    :tcol       => col  => set_color!,
    # grid mode
    :grid       => bool => set_grid!,
    # labels
    :hidelabels => bool => set_labels_off!,
    :showlabels => not  => set_labels_off!,
    :angle      => fl   => set_angle!,
    :format     => lc   => set_format!,
    :shift      => fl   => set_shift!,
    :dist       => fl   => set_dist!,
    )
merge!(TICKS_OPTS, LINESTYLE_OPTS)  # ticks line
merge!(TICKS_OPTS, TEXTSTYLE_OPTS)  # labels
set_properties!(t::Ticks; opts...) = set_properties!(TICKS_OPTS, t; opts...)

###############################################################
####
#### Options for DRAWINGS
####
###############################################################

const SCATTER2D_OPTS = OptsDict(
    :label  => str => set_labels!,
    )
merge!(SCATTER2D_OPTS, GLINESTYLE_OPTS)
merge!(SCATTER2D_OPTS, GMARKERSTYLE_OPTS)
set_properties!(s::Scatter2D; opts...) = set_properties!(SCATTER2D_OPTS, s; opts...)

const FILL2D_OPTS = OptsDict(
    :xmin  => fl  => set_xmin!,
    :xmax  => fl  => set_xmax!,
    :label => str => set_label!,
    )
merge!(FILL2D_OPTS, FILLSTYLE_OPTS)
set_properties!(f::Fill2D; opts...) = set_properties!(FILL2D_OPTS, f; opts...)

const HIST2D_OPTS = OptsDict(
    :bins    => posint => set_bins!,
    :scaling => lc     => set_scaling!,
    :horiz   => bool   => set_horiz!,
    :label   => str    => set_label!,
    )
merge!(HIST2D_OPTS, BARSTYLE_OPTS)
merge!(HIST2D_OPTS, FILLSTYLE_OPTS)
set_properties!(h::Hist2D; opts...) = set_properties!(HIST2D_OPTS, h; opts...)

const BAR2D_OPTS = OptsDict(
    :stacked => bool => set_stacked!,
    :horiz   => bool => set_horiz!,
    :label   => str  => set_labels!,
    )
merge!(BAR2D_OPTS, GBARSTYLE_OPTS)
set_properties!(gb::Bar2D; opts...) = set_properties!(BAR2D_OPTS, gb; opts...)

const BOXPLOT_OPTS = OptsDict(
    # -- global toggle
    :horiz     => bool  => set_horiz!,
    # -- now inside BOXSTYLES
    # box styling
    :bw        => posfl => set_bwidths!,
    :ww        => posfl => set_wwidths!,
    # how long should the whiskers be
    :wl        => posfl => set_wrlengths!,
    # what line style should be used to draw the boxes
    :bls       => I     => set_lstyles! ⊙ pick_blstyle,
    :blw       => posfl => set_lwidths! ⊙ pick_blstyle,
    :bc        => col   => set_colors!  ⊙ pick_blstyle,
    # median line
    :medls     => I     => set_lstyles! ⊙ pick_mlstyle,
    :medlw     => posfl => set_lwidths! ⊙ pick_mlstyle,
    :medc      => col   => set_colors!  ⊙ pick_mlstyle,

    # mean marker
    :mean_show => bool  => set_mshow!,
    :mm        => str   => set_markers! ⊙ pick_mmstyle,
    :mms       => posfl => set_msizes!  ⊙ pick_mmstyle,
    :mmc       => col   => set_mcols!   ⊙ pick_mmstyle,
    # outliers
    :out_show  => bool  => set_mshow!,
    :om        => str   => set_markers! ⊙ pick_omstyle,
    :oms       => posfl => set_msizes!  ⊙ pick_omstyle,
    :omc       => col   => set_mcols!   ⊙ pick_omstyle,
    )
set_properties!(bp::Boxplot; opts...) = set_properties!(BOXPLOT_OPTS, bp; opts...)

const HEATMAP_OPTS = OptsDict(
    :cmap     => col => set_cmap!,
    )
set_properties!(h::Heatmap; opts...) = set_properties!(HEATMAP_OPTS, h; opts...)

const SCATTER3D_OPTS = OptsDict(
    :label  => str => set_label!, # .
    )
merge!(SCATTER3D_OPTS, LINESTYLE_OPTS)
merge!(SCATTER3D_OPTS, MARKERSTYLE_OPTS)
set_properties!(s::Scatter3D; opts...) = set_properties!(SCATTER3D_OPTS, s; opts...)

###############################################################
####
#### Options for OBJECTS
####
###############################################################

const TEXT2D_OPTS = OptsDict(
    :pos => str => set_position!,
    )
merge!(TEXT2D_OPTS, TEXTSTYLE_OPTS)
set_properties!(t::Text2D; opts...) = set_properties!(TEXT2D_OPTS, t; opts...)

const STRAIGHTLINE2D_OPTS = OptsDict(
    )
merge!(STRAIGHTLINE2D_OPTS, LINESTYLE_OPTS)
set_properties!(t::StraightLine2D; opts...) = set_properties!(STRAIGHTLINE2D_OPTS, t; opts...)

const BOX2D_OPTS = OptsDict(
    :pos      => str   => set_position!,
    :nobox    => bool  => set_nobox!,
    :box      => not   => set_nobox!,
    :fill     => opcol => set_fill!,
    :alpha    => alpha => set_alpha!,
    )
merge!(BOX2D_OPTS, LINESTYLE_OPTS)
set_properties!(b::Box2D; opts...) = set_properties!(BOX2D_OPTS, b; opts...)

const COLORBAR_OPTS = OptsDict(
    :pixels     => posint => set_pixels!,
    :nobox      => bool   => set_nobox!,
    :box        => not    => set_nobox!,
    :pos        => str    => set_position!,
    :offset     => fl     => set_offset!,
    :size       => posfl  => set_size!,
    :ticks      => fl     => set_ticks!,
    :labels     => str    => set_labels!,
    )
merge!(COLORBAR_OPTS, TICKS_OPTS)
merge!(COLORBAR_OPTS, TEXTSTYLE_OPTS)
set_properties!(b::Colorbar; opts...) = set_properties!(COLORBAR_OPTS, b; opts...)

###############################################################
####
#### Options for AX*
####
###############################################################

const  AXES_OPTS = OptsDict(
    :size  => fl   => set_size!,
    :title => str  => set_title!,
    :off   => bool => set_off!,
    )
set_properties!(a::Axes; opts...) = set_properties!(AXES_OPTS, a; opts...)

const AXIS_OPTS = OptsDict(
    :title  => str   => set_title!,
    :base   => posfl => set_base!,
    :min    => fl    => set_min!,
    :max    => fl    => set_max!,
    :log    => bool  => set_log!,
    :lw     => posfl => set_lwidth!,
    :off    => bool  => set_off!,
    )
merge!(AXIS_OPTS, TEXTSTYLE_OPTS)
set_properties!(a::Axis; opts...) = set_properties!(AXIS_OPTS, a; opts...)

###############################################################
####
#### Options for FIGURE
####
###############################################################

const FIGURE_OPTS = OptsDict(
    :size     => posfl => set_size!,
    :tex      => bool  => set_texlabels!,
    :texscale => lc    => set_texscale!,
    :alpha    => bool  => set_transparency!,
    :preamble => str   => set_texpreamble!,
    :col      => opcol => set_color!,
    :bgcol    => opcol => set_color!,
    :bgalpha  => alpha => set_alpha!,
    )
merge!(FIGURE_OPTS, TEXTSTYLE_OPTS)
set_properties!(f::Figure; opts...) = set_properties!(FIGURE_OPTS, f; opts...)


text:
module Meta

using ..JOH

struct Reproducer
	n_top::Int
	n_bottom::Int
	n_random::Int
	use_mean::Bool
	select_extrema_every_perturb::Bool
	reproducer::Function
end

struct Trainer
	trainer::Function
end

struct Environment
	problem::JOH.Problem
    #algorithm behavior settings
    brood_size::Int
    population_size::Int

	#sub methods
    reproduction_method::Reproducer
    training_method::Trainer

    #runtime settings
    time_limit::Int
    num_threads::Int
end

struct Population
	solutions::Vector{JOH.Solution}
	environment::Environment
end

include("Reproduction.jl")
include("LocalSearch.jl")
include("Population.jl")

end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
function get_a(obj::DSSTThirdBodyContext)
    return jcall(obj, "getA", jdouble, ())
end

function get_alpha(obj::DSSTThirdBodyContext)
    return jcall(obj, "getAlpha", jdouble, ())
end

function get_ao_r3_pow(obj::DSSTThirdBodyContext)
    return jcall(obj, "getAoR3Pow", Vector{jdouble}, ())
end

function get_bb(obj::DSSTThirdBodyContext)
    return jcall(obj, "getBB", jdouble, ())
end

function get_bbb(obj::DSSTThirdBodyContext)
    return jcall(obj, "getBBB", jdouble, ())
end

function get_beta(obj::DSSTThirdBodyContext)
    return jcall(obj, "getBeta", jdouble, ())
end

function get_bo_a(obj::DSSTThirdBodyContext)
    return jcall(obj, "getBoA", jdouble, ())
end

function get_bo_a_bpo(obj::DSSTThirdBodyContext)
    return jcall(obj, "getBoABpo", jdouble, ())
end

function get_gamma(obj::DSSTThirdBodyContext)
    return jcall(obj, "getGamma", jdouble, ())
end

function get_hxxx(obj::DSSTThirdBodyContext)
    return jcall(obj, "getHXXX", jdouble, ())
end

function get_kxxx(obj::DSSTThirdBodyContext)
    return jcall(obj, "getKXXX", jdouble, ())
end

function get_m2ao_a(obj::DSSTThirdBodyContext)
    return jcall(obj, "getM2aoA", jdouble, ())
end

function get_m_co2ab(obj::DSSTThirdBodyContext)
    return jcall(obj, "getMCo2AB", jdouble, ())
end

function get_max_ar3_pow(obj::DSSTThirdBodyContext)
    return jcall(obj, "getMaxAR3Pow", jint, ())
end

function get_max_ecc_pow(obj::DSSTThirdBodyContext)
    return jcall(obj, "getMaxEccPow", jint, ())
end

function get_max_freq_f(obj::DSSTThirdBodyContext)
    return jcall(obj, "getMaxFreqF", jint, ())
end

function get_mean_motion(obj::DSSTThirdBodyContext)
    return jcall(obj, "getMeanMotion", jdouble, ())
end

function get_muo_r3(obj::DSSTThirdBodyContext)
    return jcall(obj, "getMuoR3", jdouble, ())
end

function get_oo_ab(obj::DSSTThirdBodyContext)
    return jcall(obj, "getOoAB", jdouble, ())
end

function get_qns(obj::DSSTThirdBodyContext)
    return jcall(obj, "getQns", Vector{Vector{jdouble}}, ())
end

function get_x(obj::DSSTThirdBodyContext)
    return jcall(obj, "getX", jdouble, ())
end

function getb(obj::DSSTThirdBodyContext)
    return jcall(obj, "getb", jdouble, ())
end


text:
module TestSKCrossValidator

using Test
using Random
using AutoMLPipeline
using AutoMLPipeline.SKCrossValidators
using AutoMLPipeline.DecisionTreeLearners
using AutoMLPipeline.Pipelines
using AutoMLPipeline.Utils

function test_skcrossvalidator()
  Random.seed!(123)
  data=getiris()
  X=data[:,1:4]
  Y=data[:,5] |> Vector{String}
  ppl1 = LinearPipeline(Dict(:machines=>[RandomForest()]))
  @test crossvalidate(ppl1,X,Y,"accuracy_score").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"balanced_accuracy_score").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"cohen_kappa_score").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"jaccard_score","weighted").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"matthews_corrcoef").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"hamming_loss").mean < 0.1
  @test crossvalidate(ppl1,X,Y,"zero_one_loss").mean < 0.1
  @test crossvalidate(ppl1,X,Y,"f1_score","weighted").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"precision_score","weighted").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"recall_score","weighted").mean > 0.80
end
@testset "CrossValidator" begin
  test_skcrossvalidator()
end

end


text:
"""Calculte the Continuous Ranked Probability Score (CRPS) of an ensemble array wrt
a true observation o. In-place version that will sort the array."""
function CRPS!(x::Array{T,1},o::T) where T
    sort!(x)
    n = length(x)
    dp = 1/n                                # probability increment per xi in x

    # find index m to split integration over x into heaviside=0 or 1
    m = findfirst(xi -> xi >= o,x)

    # integration from o to x if o is outside of the range of x
    s = m==1 ? x[1]-o : 0.0

    # integrate over the heaviside=0 part
    for i in 1:(m == nothing ? n-1 : m-1)
        s += (i*dp)^2*(x[i+1]-x[i])
    end

    # integrate over the heaviside=1 part
    for i in (m == nothing ? n : m):n-1
        s += (i*dp-1)^2*(x[i+1]-x[i])
    end

    # integration from x to o if o is outside of the range of x
    s += (m == nothing) ? o-x[end] : 0.0

    return s
end

"""Calculte the Continuous Ranked Probability Score (CRPS) of an ensemble array wrt
a true observation o."""
CRPS(x::Array{T,1},o::T) where T = CRPS!(copy(x),o)
CRPS(x::Array{T,1},o::Real) where T = CRPS!(copy(x),T(o))


text:
"""
Module Keyboard

Provides functions for controlling the keyboard by voice.

# Functions

###### Typing
- [`type`](@ref)

###### Special keys
- [`page`](@ref)

To see a description of a function type `?<functionname>`.
"""
module Keyboard

using PyCall
import ..JustSayIt: @voiceargs, pyimport_pip, controller, set_controller, TYPE_MODEL_NAME, DIGITS_ENGLISH, tic, toc, is_next, are_next, all_consumed, was_partial_recognition, InsecureRecognitionException, reset_all, do_delayed_resets

## PYTHON MODULES

const Pynput = PyNULL()

function __init__()
    ENV["PYTHON"] = ""                                              # FORCE PyCall to use Conda.jl
    copy!(Pynput, pyimport_pip("pynput"))
    set_controller("keyboard", Pynput.keyboard.Controller())
end

## CONSTANTS

const TYPE_MEMORY_ALLOC_GRANULARITY = 32
const TYPE_KEYWORDS_ENGLISH = Dict("end"=>"terminus", "undo"=>"undo", "redo"=>"redo", "uppercase"=>"uppercase", "lowercase"=>"lowercase", "digits"=>"digits", "point"=>"point", "comma"=>"comma", "colon"=>"colon", "semicolon"=>"semicolon", "exclamation"=>"exclamation", "interrogation"=>"interrogation", "paragraph" => "paragraph")

## FUNCTIONS

@doc """
    type `text` | `words` | `digits`

Type in one of the following modes:
- `text`
- `words`
- `digits`

Each of the modes supports a set of keywords which can trigger some immediate action or modify the handling of subsequent speech input. It is important to note that the speech is analysed in word groups which are naturally delimited by longer silences; and keywords are only considered as such if their word group does not contain anything else then keywords. This allows to determine whether a word that is recognised as a keyword should trigger some action or be typed instead.

# Keywords
- "terminus": end typing.
- "undo": undo typing of last word group.
- "redo": redo typing of last word group.
- "uppercase": type the first word of the next word group uppercase (automatic in `text` mode after '.', '!' and '?').
- "lowercase": type the first word of the next word group lowercase (default).
- "digits": interpret the next word group as digits (interprets in addition 'dot' as '.').
- "point": type '.'.
- "comma": type ','.
- "colon": type ':'.
- "semicolon": type ';'.
- "exclamation": type '!'.
- "interrogation": type '?'.
- "paragraph": start a new paragraph.

# Modes

#### `text`
Type any kind of text supported by the keywords, including digits and punctuation marks. All keywords are supported.

#### `words`
Type words only (upper and lower case). Supported keywords are:
- "terminus"
- "undo"
- "redo"
- "uppercase"
- "lowercase"

#### `digits`
Type digits only (including '.'). Supported keywords are:
- "terminus"
- "undo"
- "redo"
- "digits"
"""
type
@enum TokenGroupKind undefined_kind keyword_kind word_kind digit_kind punctuation_kind
@enum TypeMode text words digits
@voiceargs (mode=>(valid_input_auto=true)) function type(mode::TypeMode)
    @info "Typing $(string(mode))..."
    keyboard         = controller("keyboard")
    type_memory      = Vector{String}()
    tokengroup_str   = ""
    tokengroup_kind  = undefined_kind
    is_typing        = true
    was_keyword      = false
    is_new_group     = false
    is_uppercase     = (mode == text) ? true : false
    punctuation      = Vector{String}()
    undo_count       = 0
    ig               = 0  # group index
    it               = 0  # token index
    nb_keyword_chars = 0
    if     (mode == text)   type_keywords = [values(TYPE_KEYWORDS_ENGLISH)...]
    elseif (mode == words)  type_keywords = [TYPE_KEYWORDS_ENGLISH["end"], TYPE_KEYWORDS_ENGLISH["undo"], TYPE_KEYWORDS_ENGLISH["redo"], TYPE_KEYWORDS_ENGLISH["uppercase"], TYPE_KEYWORDS_ENGLISH["lowercase"]]
    elseif (mode == digits) type_keywords = [TYPE_KEYWORDS_ENGLISH["end"], TYPE_KEYWORDS_ENGLISH["undo"], TYPE_KEYWORDS_ENGLISH["redo"], TYPE_KEYWORDS_ENGLISH["digits"]]
    end
    while is_typing
        is_new_group = (all_consumed() && !was_partial_recognition())
        if ig == 0
            ig = 1
        elseif is_new_group && !was_keyword
            if ig > length(type_memory)
                resize!(type_memory, length(type_memory) + TYPE_MEMORY_ALLOC_GRANULARITY)
            end
            type_memory[ig] = tokengroup_str
            tokengroup_str = ""
            ig += 1
            it = 0
        end
        if is_new_group && (tokengroup_kind == undefined_kind)
            reset_all(; hard=true, exclude_active=true)      # NOTE: a reset is required to avoid that the dynamic recognizers generated in is_next and are_next include the previous recognition as desired for normal command recognition.
            if is_next(type_keywords; use_max_accuracy=false, ignore_unknown=false)
                keywords = String[]
                try
                    are_keywords, keywords = are_next(type_keywords; consume_if_match=true, ignore_unknown=false)
                catch e
                    if isa(e, InsecureRecognitionException)
                        are_keywords = false
                    else
                        rethrow(e)
                    end
                end
            else
                are_keywords = false
            end
            if are_keywords
                tokengroup_kind = keyword_kind
                reset_all(; hard=true)         # Force hard resets for performance reasons.
            else
                do_delayed_resets(; hard=true) # Force delayed resets to be done hard for performance reasons (audio must imperatively not be reset if not consumed!)
            end
        end
        if tokengroup_kind == undefined_kind
            if     (mode == text)   tokengroup_kind = word_kind
            elseif (mode == words)  tokengroup_kind = word_kind
            elseif (mode == digits) tokengroup_kind = digit_kind
            end
        end
        if tokengroup_kind == keyword_kind
            @debug "Keywords found:" keywords
            tokengroup_kind = undefined_kind
            for i = 1:length(keywords)
                keyword = keywords[i]
                if tokengroup_kind != undefined_kind
                    @info "ABORT of keyword interpretation: keyword '$(keywords[i-1])' was followed by another keyword ('$keyword'). However, keyword '$(keywords[i-1])' defines the kind of the next word group and, therefore, no other keyword can follow it."
                end
                if keyword == TYPE_KEYWORDS_ENGLISH["end"]
                    is_typing = false
                elseif keyword == TYPE_KEYWORDS_ENGLISH["undo"]
                    if ig > 1
                        ig -= 1
                        @info "Undo typing of last word group..."
                        type_backspace(;count=length(type_memory[ig]))
                        undo_count += 1
                        if (ig == 0) || (type_memory[ig-1] in [".", "!", "?"])
                            is_uppercase = true
                        end
                    else
                        @info "Nothing to undo."
                    end
                elseif keyword == TYPE_KEYWORDS_ENGLISH["redo"]
                    if undo_count > 0 && ig <= length(type_memory)
                        @info "Redo typing of last word group..."
                        if nb_keyword_chars > 0
                            type_backspace(;count=nb_keyword_chars) # NOTE: the removal of keyword signs must be done after the call to obtain the next token, in order to have it visible until it is spoken.
                            nb_keyword_chars = 0
                            is_uppercase = false
                        end
                        keyboard.type(type_memory[ig])
                        undo_count -= 1
                        ig += 1
                    else
                        @info "Nothing to redo."
                    end
                elseif keyword == TYPE_KEYWORDS_ENGLISH["uppercase"]
                    is_uppercase = true
                elseif keyword == TYPE_KEYWORDS_ENGLISH["lowercase"]
                    is_uppercase = false
                elseif keyword == TYPE_KEYWORDS_ENGLISH["digits"]
                    tokengroup_kind = digit_kind
                elseif keyword == TYPE_KEYWORDS_ENGLISH["point"]
                    push!(punctuation, ".")
                    tokengroup_kind = punctuation_kind
                    is_uppercase = true
                elseif keyword == TYPE_KEYWORDS_ENGLISH["comma"]
                    push!(punctuation, ",")
                    tokengroup_kind = punctuation_kind
                elseif keyword == TYPE_KEYWORDS_ENGLISH["colon"]
                    push!(punctuation, ":")
                    tokengroup_kind = punctuation_kind
                elseif keyword == TYPE_KEYWORDS_ENGLISH["semicolon"]
                    push!(punctuation, ";")
                    tokengroup_kind = punctuation_kind
                elseif keyword == TYPE_KEYWORDS_ENGLISH["exclamation"]
                    push!(punctuation, "!")
                    tokengroup_kind = punctuation_kind
                    is_uppercase = true
                elseif keyword == TYPE_KEYWORDS_ENGLISH["interrogation"]
                    push!(punctuation, "?")
                    tokengroup_kind = punctuation_kind
                    is_uppercase = true
                elseif keyword == TYPE_KEYWORDS_ENGLISH["paragraph"]
                    push!(punctuation, "\n")
                    tokengroup_kind = punctuation_kind
                else
                    @info "Unkown keyword." #NOTE: this should never occur as the are_next should only match known keywords.
                end
            end
            was_keyword = true
            keyword_sign = ""
            if tokengroup_kind == digit_kind
                keyword_sign = " [$(TYPE_KEYWORDS_ENGLISH["digits"])]"
            elseif is_uppercase == true
                keyword_sign = " [$(TYPE_KEYWORDS_ENGLISH["uppercase"])]"
            end
            if keyword_sign != ""
                keyboard.type(keyword_sign)
                nb_keyword_chars += length(keyword_sign)
            end
        else
            if     (tokengroup_kind == digit_kind) token = next_digit()
            elseif (tokengroup_kind == word_kind)  token = next_word()
            end
            if nb_keyword_chars > 0
                type_backspace(;count=nb_keyword_chars) # NOTE: the removal of keyword signs must be done after the call to obtain the next token, in order to have it visible until it is spoken.
                nb_keyword_chars = 0
            end
            if (tokengroup_kind == word_kind)
                if (is_uppercase) token = uppercasefirst(token) end
                if (ig == 1 && it == 0) token_str = token
                else                    token_str = " " * token
                end
            elseif (tokengroup_kind == digit_kind)
                if (ig > 1 && it == 0 && mode != digits) token_str = " " * token
                else                   token_str = token
                end
            elseif (tokengroup_kind == punctuation_kind)
                token_str = join(punctuation, "")
                punctuation = Vector{String}()
            elseif (tokengroup_kind == word_kind)
                @info "Unkown token group." #NOTE: this should never occur.
            end
            keyboard.type(token_str)
            tokengroup_str *= token_str
            it += 1
            undo_count = 0
            was_keyword = false
            if (tokengroup_kind != punctuation_kind) is_uppercase = false end
            if all_consumed() tokengroup_kind = undefined_kind end
            sleep(0.05)
        end
    end
    @info "...stopped typing $(string(mode))."
end

interpret_digits(input::AbstractString) = (return DIGITS_ENGLISH[input])

@doc "Get next digit from speech."
next_digit
@voiceargs digit=>(valid_input=[keys(DIGITS_ENGLISH)...], use_max_accuracy=true, interpret_function=interpret_digits) next_digit(digit::String) = (return digit)

@doc "Get next word from speech."
next_word
@voiceargs word=>(model=TYPE_MODEL_NAME, use_max_accuracy=true) next_word(word::String) = (return word)

function type_backspace(; count::Integer=1)
    keyboard  = controller("keyboard")
    backspace = Pynput.keyboard.Key.backspace
    for i = 1:count
        keyboard.press(backspace); keyboard.release(backspace)
    end
end

end # module Keyboard