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

struct UsersApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_usersConversations(_api::UsersApi; user=nothing, types=nothing, exclude_archived=nothing, limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/users.conversations", ["slackAuth"])
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_param(_ctx.query, "types", types)  # type String
    Swagger.set_param(_ctx.query, "exclude_archived", exclude_archived)  # type Bool
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

List conversations the calling user may access.
Param: user::String
Param: types::String
Param: exclude_archived::Bool
Param: limit::Int32
Param: cursor::String
Return: Dict{String, Any}
"""
function usersConversations(_api::UsersApi; user=nothing, types=nothing, exclude_archived=nothing, limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersConversations(_api; user=user, types=types, exclude_archived=exclude_archived, limit=limit, cursor=cursor, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersConversations(_api::UsersApi, response_stream::Channel; user=nothing, types=nothing, exclude_archived=nothing, limit=nothing, cursor=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersConversations(_api; user=user, types=types, exclude_archived=exclude_archived, limit=limit, cursor=cursor, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersDeletePhoto(_api::UsersApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/users.deletePhoto", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Delete the user profile photo
Return: Dict{String, Any}
"""
function usersDeletePhoto(_api::UsersApi; _mediaType=nothing)
    _ctx = _swaggerinternal_usersDeletePhoto(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersDeletePhoto(_api::UsersApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_usersDeletePhoto(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersGetPresence(_api::UsersApi; user=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/users.getPresence", ["slackAuth"])
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Gets user presence information.
Param: user::String
Return: Dict{String, Any}
"""
function usersGetPresence(_api::UsersApi; user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersGetPresence(_api; user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersGetPresence(_api::UsersApi, response_stream::Channel; user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersGetPresence(_api; user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersIdentity(_api::UsersApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Any, "/users.identity", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Get a user's identity.
Return: Any
"""
function usersIdentity(_api::UsersApi; _mediaType=nothing)
    _ctx = _swaggerinternal_usersIdentity(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersIdentity(_api::UsersApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_usersIdentity(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersInfo(_api::UsersApi; include_locale=nothing, user=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/users.info", ["slackAuth"])
    Swagger.set_param(_ctx.query, "include_locale", include_locale)  # type Bool
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Gets information about a user.
Param: include_locale::Bool
Param: user::String
Return: Dict{String, Any}
"""
function usersInfo(_api::UsersApi; include_locale=nothing, user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersInfo(_api; include_locale=include_locale, user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersInfo(_api::UsersApi, response_stream::Channel; include_locale=nothing, user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersInfo(_api; include_locale=include_locale, user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersList(_api::UsersApi; limit=nothing, cursor=nothing, include_locale=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/users.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "include_locale", include_locale)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Lists all users in a Slack team.
Param: limit::Int32
Param: cursor::String
Param: include_locale::Bool
Return: Dict{String, Any}
"""
function usersList(_api::UsersApi; limit=nothing, cursor=nothing, include_locale=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersList(_api; limit=limit, cursor=cursor, include_locale=include_locale, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersList(_api::UsersApi, response_stream::Channel; limit=nothing, cursor=nothing, include_locale=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersList(_api; limit=limit, cursor=cursor, include_locale=include_locale, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersLookupByEmail(_api::UsersApi, email::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/users.lookupByEmail", ["slackAuth"])
    Swagger.set_param(_ctx.query, "email", email)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Find a user with an email address.
Param: email::String (required)
Return: Dict{String, Any}
"""
function usersLookupByEmail(_api::UsersApi, email::String; _mediaType=nothing)
    _ctx = _swaggerinternal_usersLookupByEmail(_api, email; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersLookupByEmail(_api::UsersApi, response_stream::Channel, email::String; _mediaType=nothing)
    _ctx = _swaggerinternal_usersLookupByEmail(_api, email; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersProfileGet(_api::UsersApi; include_labels=nothing, user=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/users.profile.get", ["slackAuth"])
    Swagger.set_param(_ctx.query, "include_labels", include_labels)  # type Bool
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Retrieves a user's profile information.
Param: include_labels::Bool
Param: user::String
Return: Dict{String, Any}
"""
function usersProfileGet(_api::UsersApi; include_labels=nothing, user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersProfileGet(_api; include_labels=include_labels, user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersProfileGet(_api::UsersApi, response_stream::Channel; include_labels=nothing, user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersProfileGet(_api; include_labels=include_labels, user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersProfileSet(_api::UsersApi; name=nothing, profile=nothing, user=nothing, value=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/users.profile.set", ["slackAuth"])
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_param(_ctx.form, "profile", profile)  # type String
    Swagger.set_param(_ctx.form, "user", user)  # type String
    Swagger.set_param(_ctx.form, "value", value)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set the profile information for a user.
Param: name::String
Param: profile::String
Param: user::String
Param: value::String
Return: Dict{String, Any}
"""
function usersProfileSet(_api::UsersApi; name=nothing, profile=nothing, user=nothing, value=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersProfileSet(_api; name=name, profile=profile, user=user, value=value, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersProfileSet(_api::UsersApi, response_stream::Channel; name=nothing, profile=nothing, user=nothing, value=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersProfileSet(_api; name=name, profile=profile, user=user, value=value, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersSetActive(_api::UsersApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/users.setActive", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Marked a user as active. Deprecated and non-functional.
Return: Dict{String, Any}
"""
function usersSetActive(_api::UsersApi; _mediaType=nothing)
    _ctx = _swaggerinternal_usersSetActive(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersSetActive(_api::UsersApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_usersSetActive(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersSetPhoto(_api::UsersApi; crop_w=nothing, crop_x=nothing, crop_y=nothing, image=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/users.setPhoto", ["slackAuth"])
    Swagger.set_param(_ctx.form, "crop_w", crop_w)  # type String
    Swagger.set_param(_ctx.form, "crop_x", crop_x)  # type String
    Swagger.set_param(_ctx.form, "crop_y", crop_y)  # type String
    Swagger.set_param(_ctx.form, "image", image)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Set the user profile photo
Param: crop_w::String
Param: crop_x::String
Param: crop_y::String
Param: image::String
Return: Dict{String, Any}
"""
function usersSetPhoto(_api::UsersApi; crop_w=nothing, crop_x=nothing, crop_y=nothing, image=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersSetPhoto(_api; crop_w=crop_w, crop_x=crop_x, crop_y=crop_y, image=image, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersSetPhoto(_api::UsersApi, response_stream::Channel; crop_w=nothing, crop_x=nothing, crop_y=nothing, image=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_usersSetPhoto(_api; crop_w=crop_w, crop_x=crop_x, crop_y=crop_y, image=image, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_usersSetPresence(_api::UsersApi, presence::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/users.setPresence", ["slackAuth"])
    Swagger.set_param(_ctx.form, "presence", presence)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Manually sets user presence.
Param: presence::String (required)
Return: Dict{String, Any}
"""
function usersSetPresence(_api::UsersApi, presence::String; _mediaType=nothing)
    _ctx = _swaggerinternal_usersSetPresence(_api, presence; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function usersSetPresence(_api::UsersApi, response_stream::Channel, presence::String; _mediaType=nothing)
    _ctx = _swaggerinternal_usersSetPresence(_api, presence; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export usersConversations, usersDeletePhoto, usersGetPresence, usersIdentity, usersInfo, usersList, usersLookupByEmail, usersProfileGet, usersProfileSet, usersSetActive, usersSetPhoto, usersSetPresence


text:
struct CircularVector{T}
    v::Vector{T}
    l::Int
    CircularVector{T}(l::Integer) where T = new(zeros(T, l), l)
end
function min_abs_diff(v::CircularVector{T}) where T
    val = Inf
    for i in 1:Base.length(v)
        val = min(val, abs(v[i] - v[i-1]))
    end
    return val
end
function max_abs_diff(v::CircularVector{T}) where T
    val = 0.0
    for i in 1:Base.length(v)
        val = max(val, abs(v[i] - v[i-1]))
    end
    return val
end

function Base.getindex(V::CircularVector{T}, i::Int) where T
    return V.v[mod1(i, V.l)]
end
function Base.setindex!(V::CircularVector{T}, val::T, i::Int) where T
    V.v[mod1(i, V.l)] = val
end
function Base.length(V::CircularVector{T}) where T
    return V.l
end

Base.@kwdef mutable struct Options

    # Printing options
    log_verbose::Bool = false
    log_freq::Int = 1000
    timer_verbose::Bool = false
    timer_file::Bool = false
    disable_julia_logger = true

    # time options
    time_limit::Float64 = 3600_00. #100 hours

    warn_on_limit::Bool = false
    extended_log::Bool = false
    extended_log2::Bool = false
    log_repeat_header::Bool = false

    # Default tolerances
    tol_gap::Float64 = 1e-4
    tol_feasibility::Float64 = 1e-4
    tol_feasibility_dual::Float64 = 1e-4
    tol_primal::Float64 = 1e-4
    tol_dual::Float64 = 1e-4
    tol_psd::Float64 = 1e-7
    tol_soc::Float64 = 1e-7

    check_dual_feas::Bool = false
    check_dual_feas_freq::Int = 1000

    max_obj::Float64 = 1e20
    min_iter_max_obj::Int = 10

    # infeasibility check
    min_iter_time_infeas::Int = 1000
    infeas_gap_tol::Float64 = 1e-4
    infeas_limit_gap_tol::Float64 = 1e-1
    infeas_stable_gap_tol::Float64 = 1e-4
    infeas_feasibility_tol::Float64 = 1e-4
    infeas_stable_feasibility_tol::Float64 = 1e-8

    certificate_search::Bool = true
    certificate_obj_tol::Float64 = 1e-1
    certificate_fail_tol::Float64 = 1e-8

    # Bounds on beta (dual_step / primal_step) [larger bounds may lead to numerical inaccuracy]
    min_beta::Float64 = 1e-5
    max_beta::Float64 = 1e+5
    initial_beta::Float64 = 1.

    # Adaptive primal-dual steps parameters [adapt_decay above .7 may lead to slower convergence]
    initial_adapt_level::Float64 = .9
    adapt_decay::Float64 = .8
    adapt_window::Int64 = 50

    # PDHG parameters
    convergence_window::Int = 200
    convergence_check::Int = 50
    max_iter::Int = 0
    min_iter::Int = 40
    divergence_min_update::Int = 50
    max_iter_lp::Int = 10_000_000
    max_iter_conic::Int = 1_000_000
    max_iter_local::Int = 0 #ignores user setting

    advanced_initialization::Bool = true

    # Linesearch parameters
    line_search_flag::Bool = true
    max_linsearch_steps::Int = 5000
    delta::Float64 = .9999
    initial_theta::Float64 = 1.
    linsearch_decay::Float64 = .75

    # Spectral decomposition parameters
    full_eig_decomp::Bool = false
    max_target_rank_krylov_eigs::Int = 16
    min_size_krylov_eigs::Int = 100
    warm_start_eig::Bool = true
    rank_increment::Int = 1 # 0=multiply, 1 = add
    rank_increment_factor::Int = 1 # 0 multiply, 1 = add

    # eigsolver selection
    #=
        1: Arpack [dsaupd] (tipically non-deterministic)
        2: KrylovKit [eigsolve/Lanczos] (DEFAULT)
    =#
    eigsolver::Int = 2
    eigsolver_min_lanczos::Int = 25
    eigsolver_resid_seed::Int = 1234

    # Arpack
    # note that Arpack is Non-deterministic
    # (https://github.com/mariohsouto/ProxSDP.jl/issues/69)
    arpack_tol::Float64 = 1e-10
    #=
        0: arpack random [usually faster - NON-DETERMINISTIC - slightly]
        1: all ones [???]
        2: julia random uniform (eigsolver_resid_seed) [medium for DETERMINISTIC]
        3: julia normalized random normal (eigsolver_resid_seed) [best for DETERMINISTIC]
    =#
    arpack_resid_init::Int = 3
    arpack_reset_resid::Bool = true # true for determinism
    # larger is more stable to converge and more deterministic
    arpack_max_iter::Int = 10_000
    # see remark for of dsaupd

    # KrylovKit
    krylovkit_reset_resid::Bool = false
    krylovkit_resid_init::Int = 3
    krylovkit_tol::Float64 = 1e-12
    krylovkit_max_iter::Int = 100
    krylovkit_eager::Bool = false
    krylovkit_verbose::Int = 0

    # Reduce rank [warning: heuristics]
    reduce_rank::Bool = false
    rank_slack::Int = 3

    full_eig_freq::Int = 10000000
    full_eig_len::Int = 0

    # equilibration parameters
    equilibration::Bool = false
    equilibration_iters::Int = 1000
    equilibration_lb::Float64 = -10.0
    equilibration_ub::Float64 = +10.0
    equilibration_limit::Float64 = 0.9
    equilibration_force::Bool = false

    # spectral norm [using exact norm via svds may result in nondeterministic behavior]
    approx_norm::Bool = true
end

mutable struct AffineSets
    n::Int  # Size of primal variables
    p::Int  # Number of linear equalities
    m::Int  # Number of linear inequalities
    extra::Int  # Number of adition linear equalities (for disjoint cones)
    A::SparseMatrixCSC{Float64,Int64}
    G::SparseMatrixCSC{Float64,Int64}
    b::Vector{Float64}
    h::Vector{Float64}
    c::Vector{Float64}
end

mutable struct SDPSet
    vec_i::Vector{Int}
    mat_i::Vector{Int}
    tri_len::Int
    sq_len::Int
    sq_side::Int
end

mutable struct SOCSet
    idx::Vector{Int}
    len::Int
end

mutable struct ConicSets
    sdpcone::Vector{SDPSet}
    socone::Vector{SOCSet}
end

mutable struct CPResult
    status::Int
    status_string::String
    primal::Vector{Float64}
    dual_cone::Vector{Float64}
    dual_eq::Vector{Float64}
    dual_in::Vector{Float64}
    slack_eq::Vector{Float64}
    slack_in::Vector{Float64}
    primal_residual::Float64
    dual_residual::Float64
    objval::Float64
    dual_objval::Float64
    gap::Float64
    time::Float64
    final_rank::Int
    primal_feasible_user_tol::Bool
    dual_feasible_user_tol::Bool
    certificate_found::Bool
end

mutable struct PrimalDual
    x::Vector{Float64}
    x_old::Vector{Float64}
    y::Vector{Float64}
    y_old::Vector{Float64}

    PrimalDual(aff::AffineSets) = new(
        zeros(aff.n), zeros(aff.n), zeros(aff.m+aff.p), zeros(aff.m+aff.p)
    )
end

mutable struct WarmStart
    x::Vector{Float64}
    y_eq::Vector{Float64}
    y_in::Vector{Float64}
end

mutable struct Residuals
    dual_gap::CircularVector{Float64}
    prim_obj::CircularVector{Float64}
    dual_obj::CircularVector{Float64}
    equa_feasibility::Float64
    ineq_feasibility::Float64
    feasibility::CircularVector{Float64}
    primal_residual::CircularVector{Float64}
    dual_residual::CircularVector{Float64}
    comb_residual::CircularVector{Float64}

    Residuals(window::Int) = new(
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window),
        .0,
        .0,
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window)
    )
end

const ViewVector = SubArray#{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int}}, true}
const ViewScalar = SubArray#{Float64, 1, Vector{Float64}, Tuple{Int}, true}

mutable struct AuxiliaryData
    m::Vector{Symmetric{Float64,Matrix{Float64}}}

    Mty::Vector{Float64}
    Mty_old::Vector{Float64}

    Mx::Vector{Float64}
    Mx_old::Vector{Float64}

    y_half::Vector{Float64}
    y_temp::Vector{Float64}

    soc_v::Vector{ViewVector}
    soc_s::Vector{ViewScalar}

    function AuxiliaryData(aff::AffineSets, cones::ConicSets)
        new(
            [Symmetric(zeros(sdp.sq_side, sdp.sq_side), :U) for sdp in cones.sdpcone],
            zeros(aff.n), zeros(aff.n),
            zeros(aff.p+aff.m), zeros(aff.p+aff.m),
            zeros(aff.p+aff.m), zeros(aff.p+aff.m),
            ViewVector[], ViewScalar[]
        )
    end
end

mutable struct Matrices
    M::SparseMatrixCSC{Float64,Int64}
    Mt::SparseMatrixCSC{Float64,Int64}
    c::Vector{Float64}

    Matrices(M, Mt, c) = new(M, Mt, c)
end

mutable struct Params
    current_rank::Vector{Int}
    target_rank::Vector{Int}
    rank_update::Int
    update_cont::Int
    min_eig::Vector{Float64}
    iter::Int
    stop_reason::Int
    stop_reason_string::String
    iteration::Int
    primal_step::Float64
    primal_step_old::Float64
    dual_step::Float64
    theta::Float64
    beta::Float64
    adapt_level::Float64
    window::Int
    time0::Float64
    norm_c::Float64
    norm_b::Float64
    norm_h::Float64
    sqrt2::Float64

    dual_feasibility::Float64
    dual_feasibility_check::Bool

    certificate_search::Bool
    certificate_search_min_iter::Int
    certificate_found::Bool

    # solution backup

    Params() = new()
end

ext: jl


text:
function duration(vs,fs,d) #That function cut samples length of any sound according to their duration.
    if(d<=size(vs)[1]/fs-1/fs)
    k=int64(d*fs+1)
    Vs=zeros(k,size(vs)[2])
    for j=1:size(vs)[2]
        for i=1:size(vs)[1]
            if(i<k+1)
        Vs[i,j]=vs[i,j]
            end
        end
    end
    return Vs
    else
        error("Tiempo excedido")
    end
end
#Inicio del programa

dur=5; # duration of sound signal
using WAV #That package let me read any sound signal and convert that signal in samples
using PyPlot #Plot Pakage
using DSP #function for Digital Signal Processing
close("all")

archivo1="C:\\JULIA\\Grab6.wav"
vs1, fs1 = wavread(archivo1)
Vs1=duration(vs1,fs1,dur)
Vs1
length((Vs1))
archivo2="C:\\JULIA\\Grab7.wav"
vs2, fs2 = wavread(archivo2)
Vs2=duration(vs2,fs2,dur)
Vs2
archivo3="C:\\JULIA\\Grab8.wav"
vs3, fs3 = wavread(archivo3)
Vs3=duration(vs3,fs3,dur)
Vs3

dt=1/fs1;
#te=[dt:dt:dt*(size(Vs1)[1])]
t=linspace(dt,dt*size(Vs1)[1],size(Vs1)[1])
plot(t,Vs1[:,1])
##plot(te,Vs1[:,1])
ruido1=0.006/20*rand(size(t))
#Señales con ruido
x1=Vs1[:,1]+ruido1
x2=Vs2[:,1]+ruido1
x3=Vs3[:,1]+ruido1

f=linspace(-1*fs1/2,1*fs1/2,length(t))
figure(1)
espectro1=abs(fftshift(fft(x1)/length(t)))
espectro2=abs(fftshift(fft(x2)/length(t)))
espectro3=abs(fftshift(fft(x3)/length(t)))
subplot(3,2,1)
plot(f,espectro1)
title("Primera señal")
subplot(3,2,3)
plot(f,espectro2)
title("Segunda señal")
subplot(3,2,5)
plot(f,espectro3)
title("Tercera señal")
#filtrado de señal con ruido

#Filtracion del ruido

    #Primera señal
response1=Lowpass(3000,fs=fs1)
orden=Butterworth(10);
filtrada1=filt(digitalfilter(response1,orden),x1)
FILTRADA1=abs(fftshift(fft(filtrada1)/length(t)))
subplot(3,2,2)
plot(f,FILTRADA1)
title("Primera señal filtrada 3Khz")

    #Segunda señal
response2=Lowpass(3000,fs=fs2)
orden=Butterworth(10);
filtrada2=filt(digitalfilter(response2,orden),x2)
FILTRADA2=abs(fftshift(fft(filtrada2)/length(t)))
subplot(3,2,4)
plot(f,FILTRADA2)
title("Segunda señal filtrada 3Khz")

    #Tercera señal
response3=Lowpass(4000,fs=fs3)
orden=Butterworth(10);
filtrada3=filt(digitalfilter(response3,orden),x3)
FILTRADA3=abs(fftshift(fft(filtrada3)/length(t)))
subplot(3,2,6)
plot(f,FILTRADA3)
title("Tercera señal filtrada 4Khz")

    #ModulacionAM Doble Banda lateral
figure(2)
title("ModulacionAM")
dbl2=filtrada2.*cos(2*pi*3500*t)
DBL2=abs(fftshift(fft(dbl2)/length(t)))
subplot(2,2,1)
plot(f,DBL2)
title("Modulacion 3.5Khz segunda señal (DBL)")
response22=Bandpass(3500,6500,fs=fs2)
blu2=filt(digitalfilter(response22,orden),dbl2)
BLU2=abs(fftshift(fft(blu2)/length(t)))
subplot(2,2,2)
plot(f,BLU2)
title("BLU segunda señal")
dbl3=filtrada3.*cos(2*pi*7000*t)
DBL3=abs(fftshift(fft(dbl3)/length(t)))
subplot(2,2,3)
plot(f,DBL3)
title("Modulacion 7Khz tercera señal (DBL)")
response33=Bandpass(7000,11000,fs=fs3)
blu3=filt(digitalfilter(response33,orden),dbl3)
BLU3=abs(fftshift(fft(blu3)/length(t)))
subplot(2,2,4)
plot(f,BLU3)
title("BLU tercera señal")
figure(3)
plot(f,FILTRADA1)
plot(f,BLU2)
plot(f,BLU3)
title("Señal FDM")
    #Señal FDM con ruido
suma=filtrada1+blu2+blu3+ruido1
FDM=abs(fftshift(fft(suma)/length(t)))
figure(4)
subplot(3,1,1)
plot(f,FDM)
title("Señal FDM con ruido")
    #Filtracion de señal FDM
response_FDM=Lowpass(11000,fs=fs1)
orden_FDM=Butterworth(15)
fil_FDM=filt(digitalfilter(response_FDM,orden_FDM),suma)
FIL_FDM=abs(fftshift(fft(fil_FDM)/length(t)))
subplot(3,1,2)
plot(f,FIL_FDM)
title("Señal FDM filtrada")
    #ModulacionFDM
mod_FDM=fil_FDM.*cos(2*pi*11100*t)
MOD_FDM=abs(fftshift(fft(mod_FDM)/length(t)))
subplot(3,1,3)
plot(f,MOD_FDM)
title("FDM modulada a 11.2Khz")

max_stars_repo_name: OscarNeciosup/AmplitudeModulation

max_issues_repo_name: OscarNeciosup/AmplitudeModulation

max_forks_repo_name: OscarNeciosup/AmplitudeModulation


text:
macro pywrap(ex)
    if @capture(ex, struct T_ fields__ end)
        def = MacroTools.splitstructdef(ex)
        pushfirst!(def[:fields], (:o, :(PyCall.PyObject)))
        fnames = Tuple(first(f) for f in def[:fields])
        return quote
            $(esc(MacroTools.combinestructdef(def)))
            function $Base.getproperty(obj::$(esc(def[:name])), sym::Symbol)
                if sym in ($(map(QuoteNode, fnames)...), )
                    getfield(obj, sym)
                else
                    $Base.getproperty(getfield(obj, :o), sym)
                end
            end
            function $Base.setproperty!(obj::$(esc(def[:name])), sym::Symbol, val)
                if sym in ($(map(QuoteNode, fnames)...), )
                    setfield!(obj, sym, val)
                else
                    $Base.setproperty!(getfield(obj, :o), sym, val)
                end
            end
        end
    else
        error("Expected a struct definition, got: $ex")
    end
end


text:
export TransferOperatorApproximation
import Printf: @sprintf

"""
    TransferOperatorApproximation(g, M, params)

An approximation to the transfer operator over some partition.

The transfer operator generator `g` contains the input data,
initialization steps, partition specification, and the transition
probability estimation method.
Along with the additional `params` (a named tuple), this information
is used to compute the transfer matrix `M`.

`M`  is a row-stochastic (and partition-dependent)
Markov matrix containing the approximated transition probabilities
between partition elements. If the input data were an orbit of
some dynamical system, then every partition element corresponds to
a unique (discretized) state of the system. Hence, the element ``M_{ij}``
gives the probability of jumping between state ``i`` to state ``j`` in
any particular time step.
"""
struct TransferOperatorApproximation{G <: TransferOperatorGenerator, T <: Real, P <: NamedTuple}
    g::G
    M::AbstractArray{T, 2}
    params::P
end

function Base.summary(to::T) where T<:TransferOperatorApproximation
    method = to.g.method
    data = to.g.pts
    s = size(to.M)
    l = length(to.M)
    percent_nonzero = @sprintf("%.4f", count(to.M .> 0.0)/length(to.M) * 100)

    sM = "  Transfer matrix: size $s with $percent_nonzero% nonzero entries"
    sdata = "  Input data: $data"
    smethod = "  Partition and estimation method: $method"
    return "TransferOperatorApproximation\n$(sM)\n$(sdata)\n$(smethod)"
end

Base.show(io::IO, to::T) where {T<:TransferOperatorApproximation} = println(io, summary(to))


text:
# TRR molecular trajectory file format

module TRR

using ..Dorothy
using ..Dorothy.Geometry
using ..Dorothy.PBC
using ..Dorothy.XDR
using ..Dorothy.Utils
using Formats
using FormatStreams

export register_trr, readtrr!, writetrr

const trr_format_signature = "\0\0\a\xc9\0\0\0\r\0\0\0\fGMX_trn_file"
const trr_precision_types = [Float32, Float64]

function register_trr()
	Formats.addformat("trajectory/x-trr")
	Formats.addextension("trajectory/x-trr", ".trr")
	Formats.addsignature("trajectory/x-trr", trr_format_signature)
	Formats.addreader("trajectory/x-trr", DorothyIO())
	Formats.addwriter("trajectory/x-trr", DorothyIO())
	FormatStreams.addstreamer("trajectory/x-trr", DorothyIO())
end

FormatStreams.streamf(::DorothyIO, ::MIME"trajectory/x-trr", io::T, args...;
		kwargs...) where {T<:IO} = TRRStream{T}(io, args...; kwargs...)

Base.read(::DorothyIO, ::MIME"trajectory/x-trr", io::IO, args...; kwargs...) =
		readtrr!(io, MolecularModel(), args...; kwargs...)

Base.read!(::DorothyIO, ::MIME"trajectory/x-trr", io::IO, model::MolecularModel,
		args...; kwargs...) = readtrr!(io, model, args...; kwargs...)

Base.write(::DorothyIO, ::MIME"trajectory/x-trr", io::IO,
		model::ParticleCollection, args...; kwargs...) =
		writetrr(io, model, args...; kwargs...)

struct TRRMetadata
	irsize::Int
	esize::Int
	cellsize::Int
	virialsize::Int
	pressuresize::Int
	topsize::Int
	symsize::Int
	Rsize::Int
	Vsize::Int
	Fsize::Int
	nparticles::Int
	precision::Type
	framesize::Int

	function TRRMetadata(irsize::Integer, esize::Integer, cellsize::Integer,
			virialsize::Integer, pressuresize::Integer, topsize::Integer,
			symsize::Integer, Rsize::Integer, Vsize::Integer,
			Fsize::Integer, nparticles::Integer, precision::Type)
		framesize = 76 + 2 * sizeof(precision) + cellsize + virialsize +
				pressuresize + Rsize + Vsize + Fsize
		new(irsize, esize, cellsize, virialsize, pressuresize, topsize,
				symsize, Rsize, Vsize, Fsize, nparticles, precision,
				framesize)
	end
end

function read_trr_metadata(io::IO)
	skip(io, 24)
	irsize = read_xdr_int(io)
	esize = read_xdr_int(io)
	cellsize = read_xdr_int(io)
	virialsize = read_xdr_int(io)
	pressuresize = read_xdr_int(io)
	topsize = read_xdr_int(io)
	symsize = read_xdr_int(io)
	Rsize = read_xdr_int(io)
	Vsize = read_xdr_int(io)
	Fsize = read_xdr_int(io)
	nparticles = read_xdr_int(io)
	precision = guess_trr_precision([cellsize, virialsize, pressuresize],
			[Rsize, Vsize, Fsize], nparticles)
	TRRMetadata(irsize, esize, cellsize, virialsize, pressuresize,
			topsize, symsize, Rsize, Vsize, Fsize, nparticles,
			precision)
end

function peek_trr_metadata(io::IO)
	mark(io)
	meta = read_trr_metadata(io)
	reset(io)
	meta
end

skip_trr_metadata(io::IO) = skip(io, 68)

function write_trr_metadata(io::IO, meta::TRRMetadata)
	write(io, trr_format_signature)
	write_xdr_int(io, meta.irsize)
	write_xdr_int(io, meta.esize)
	write_xdr_int(io, meta.cellsize)
	write_xdr_int(io, meta.virialsize)
	write_xdr_int(io, meta.topsize)
	write_xdr_int(io, meta.symsize)
	write_xdr_int(io, meta.pressuresize)
	write_xdr_int(io, meta.Rsize)
	write_xdr_int(io, meta.Vsize)
	write_xdr_int(io, meta.Fsize)
	write_xdr_int(io, meta.nparticles)
end

function trr_metadata_for_model(model::MolecularModel)
	properties = Symbol[]
	for key in [:cell, :virial, :pressure]
		if haskey(model.header, key)
			push!(properties, key)
		end
	end
	for key in [:R, :V, :F]
		if haskey(model, key)
			push!(properties, key)
		end
	end
	trr_metadata_from_vals(length(model), Float32, properties)
end

function trr_metadata_from_vals(nparticles::Integer, precision::Type = Float32,
		properties::AbstractVector{Symbol} = [:cell, :R])
	@boundscheck begin
		nparticles >= 0 || error("expected positive number of particles")
		precision in trr_precision_types ||
				error("invalid TRR trajectory precision")
	end
	cellsize = (:cell in properties ? 3*3 * sizeof(precision) : 0)
	virialsize = (:virial in properties ? 3*3 * sizeof(precision) : 0)
	pressuresize = (:pressure in properties ? 3*3 * sizeof(precision) : 0)
	Rsize = (:R in properties ? 3*nparticles * sizeof(precision) : 0)
	Vsize = (:V in properties ? 3*nparticles * sizeof(precision) : 0)
	Fsize = (:F in properties ? 3*nparticles * sizeof(precision) : 0)
	TRRMetadata(0, 0, cellsize, virialsize, pressuresize, 0, 0, Rsize,
			Vsize, Fsize, nparticles, precision)
end

function guess_trr_precision(header_array_sizes::AbstractVector{<:Integer},
		body_array_sizes::AbstractVector{<:Integer}, nparticles::Integer)
	for array_size in header_array_sizes
		for T in trr_precision_types
			if array_size == 3 * 3 * sizeof(T)
				return T
			end
		end
	end
	for array_size in body_array_sizes
		for T in trr_precision_types
			if array_size == 3 * nparticles * sizeof(T)
				return T
			end
		end
	end
	error("could not determine TRR trajectory precision")
end

struct TRRBuffer{T<:Real}
	cell::Matrix{T}
	cellf::Matrix{Float64}
	virial::Matrix{T}
	virialf::Matrix{Float64}
	pressure::Matrix{T}
	pressuref::Matrix{Float64}
	R::Matrix{T}
	V::Matrix{T}
	F::Matrix{T}

	function TRRBuffer{T}(meta::TRRMetadata) where {T<:Real}
		n = meta.cellsize > 0 ? 3 : 0
		cell = Matrix{T}(undef, 3,n)
		cellf = Matrix{Float64}(undef, 3,n)
		n = meta.virialsize > 0 ? 3 : 0
		virial = Matrix{T}(undef, 3,n)
		virialf = Matrix{Float64}(undef, 3,n)
		n = meta.pressuresize > 0 ? 3 : 0
		pressure = Matrix{T}(undef, 3,n)
		pressuref = Matrix{Float64}(undef, 3,n)
		n = meta.Rsize > 0 ? meta.nparticles : 0
		R = Matrix{T}(undef, 3,n)
		n = meta.Vsize > 0 ? meta.nparticles : 0
		V = Matrix{T}(undef, 3,n)
		n = meta.Fsize > 0 ? meta.nparticles : 0
		F = Matrix{T}(undef, 3,n)
		new(cell, cellf, virial, virialf, pressure, pressuref, R, V, F)
	end
end

TRRBuffer(meta::TRRMetadata) = TRRBuffer{meta.precision}(meta)

mutable struct TRRStream{T<:IO} <: MolecularTrajectory
	io::T
	origin::Int
	pos::Int
	nframes::Int
	meta::TRRMetadata
	buffer::TRRBuffer

	function TRRStream{T}(io::T; nparticles::Integer = -1,
			precision::Type = Float32,
			properties::AbstractVector{Symbol} = [:cell, :R]) where {T<:IO}
		origin = position(io)
		if eof(io)
			nparticles >= 0 || error("expected positive number of particles")
			meta = trr_metadata_from_vals(nparticles, precision, properties)
		else
			meta = read_trr_metadata(io)
		end
		seekend(io)
		nframes = (position(io) - origin) ÷ meta.framesize
		seek(io, origin)
		new(io, origin, 0, nframes, meta, TRRBuffer(meta))
	end
end

Base.close(s::TRRStream) = close(s.io)

Base.length(s::TRRStream) = s.nframes

Base.eof(s::TRRStream) = (s.pos >= s.nframes)

Base.position(s::TRRStream) = s.pos

function Base.seek(s::TRRStream, pos::Integer)
	@boundscheck checkbounds(0:s.nframes, pos)
	seek(s.io, s.origin + s.meta.framesize * pos)
	s.pos = pos
	s
end

function Base.seekstart(s::TRRStream)
	seek(s.io, s.origin)
	s.pos = 0
	s
end

function Base.seekend(s::TRRStream)
	seekend(s.io)
	s.pos = s.nframes
	s
end

Base.read(s::TRRStream) = read!(s, MolecularModel())

function Base.read!(s::TRRStream, model::MolecularModel)
	readtrr!(s.io, model, s.meta, s.buffer)
	s.pos += 1
	model
end

function Base.write(s::TRRStream, model::ParticleCollection)
	pos0 = position(s.io)
	writetrr(s.io, model, s.meta, s.buffer)
	s.pos += 1
	s.nframes += 1
	position(s.io) - pos0
end

function Base.truncate(s::TRRStream, n::Integer)
	@boundscheck n >= 0 || error("expected positive value")
	truncate(s.io, n * s.meta.framesize)
	s.nframes = n
	if s.pos > n
		s.pos = n
	end
	s
end

function readtrr!(io::IO, model::MolecularModel,
		meta::TRRMetadata = peek_trr_metadata(io),
		buffer::TRRBuffer = TRRBuffer(meta))
	skip_trr_metadata(io)
	resize!(model, meta.nparticles)
	model.header.step = read_xdr_int(io)
	skip(io, 4)
	model.header.time = read_xdr_num(io, meta.precision)
	model.header.lambda = read_xdr_num(io, meta.precision)
	if meta.cellsize != 0
		read!(io, buffer.cell)
		decode_trr_array!(buffer.cellf, buffer.cell, 10.0)
		model.header.cell = TriclinicCell(buffer.cellf)
	end
	if meta.virialsize != 0
		read!(io, buffer.virial)
		decode_trr_array!(buffer.virialf, buffer.virial, 1.0)
		model.header.virial = buffer.virialf
	end
	if meta.pressuresize != 0
		read!(io, buffer.pressure)
		decode_trr_array!(buffer.pressuref, buffer.pressure, 1.0/1000.0)
		model.header.pressure = buffer.pressuref
	end
	if meta.Rsize != 0
		read!(io, buffer.R)
		decode_trr_array!(get!(model, :R, undef), buffer.R, 10.0)
	end
	if meta.Vsize != 0
		read!(io, buffer.V)
		decode_trr_array!(get!(model, :V, undef), buffer.V, 10.0)
	end
	if meta.Fsize != 0
		read!(io, buffer.F)
		decode_trr_array!(get!(model, :F, undef), buffer.F, 1/10.0)
	end
	model
end

function decode_trr_array!(dest::Array{Float64}, src::Array{<:Real},
		factor::Real)
	@inbounds @fastmath @simd for i in eachindex(dest)
		dest[i] = ntoh(src[i]) * factor
	end
	dest
end

function decode_trr_array!(dest::Vector{Vector3D}, src::Array{<:Real},
		factor::Real)
	@inbounds @fastmath @simd for i in eachindex(dest)
		d = (i-1) * 3
		x = ntoh(src[d+1]) * factor
		y = ntoh(src[d+2]) * factor
		z = ntoh(src[d+3]) * factor
		dest[i] = Vector3D(x,y,z)
	end
	dest
end

function writetrr(io::IO, model::ParticleCollection,
		meta::TRRMetadata = trr_metadata_for_model(model),
		buffer::TRRBuffer = TRRBuffer(meta))
	write_trr_metadata(io, meta)
	write_xdr_int(io, model.header.step)
	write_xdr_int(io, 0)
	write_xdr_num(io, model.header.time, meta.precision)
	write_xdr_num(io, model.header.lambda, meta.precision)
	if meta.cellsize != 0
		buffer.cellf .= pbcmatrix(model.header.cell)
		encode_trr_array!(buffer.cell, buffer.cellf, 10.0)
		write(io, buffer.cell)
	end
	if meta.virialsize != 0
		buffer.virialf .= model.header.virial
		encode_trr_array!(buffer.virial, buffer.virialf, 1.0)
		write(io, buffer.virial)
	end
	if meta.pressuresize != 0
		buffer.pressuref .= model.header.pressure
		encode_trr_array!(buffer.pressure, buffer.pressuref, 1.0/1000.0)
		write(io, buffer.pressure)
	end
	if meta.Rsize != 0
		encode_trr_array!(buffer.R, model.R, 10.0)
		write(io, buffer.R)
	end
	if meta.Vsize != 0
		encode_trr_array!(buffer.V, model.V, 10.0)
		write(io, buffer.V)
	end
	if meta.Fsize != 0
		encode_trr_array!(buffer.F, model.F, 1/10.0)
		write(io, buffer.F)
	end
end

function encode_trr_array!(dest::Array{T}, src::AbstractArray{Float64},
		factor::Real) where {T<:Real}
	@inbounds @fastmath @simd for i in eachindex(src)
		dest[i] = hton(T(src[i] / factor))
	end
	dest
end

function encode_trr_array!(dest::Array{T}, src::AbstractVector{Vector3D},
		factor::Real) where {T<:Real}
	@inbounds @fastmath @simd for i in eachindex(src)
		x, y, z = hton.(T.(src[i] ./ factor))
		d = (i-1) * 3
		dest[d+1] = x
		dest[d+2] = y
		dest[d+3] = z
	end
	dest
end

end # module


text:
function naiveLU( a:: Array{Float64,1}, b :: Array{Float64,1}, c :: Array{Float64,1})
  #= @assert (size(a) == size(b)) =#
  #= @assert (size(a) == size(c)) =#
  l = zeros(a)
  d = zeros(a)
  u = zeros(a)
  d[1] = a[1]
  u[1:end-1] = c[1:end-1]

  for ii = 2:(length(l))
    l[ii] = b[ii]/d[ii-1]
    d[ii] = a[ii] - l[ii]u[ii-1]
  end
  return (l,d,u)
end
export naiveLU

function ForwardThomas( l :: Array{Float64,1}, dataR :: Array{Float64,1})
  y = zeros(dataR)
  y[1] = dataR[1]
  for ii = 2:length(dataR)
    y[ii] = dataR[ii] - l[ii]*y[ii-1]
  end
  return y
end

function BackwardThomas( d :: Array{Float64, 1}, u :: Array{Float64,1}, dataR :: Array{Float64,1})
  x = zeros(dataR)
  x[end] = dataR[end]/d[end]
  for ii = (length(x)-1):-1:1
    x[ii] = (dataR[ii] - u[ii]*x[ii+1])/d[ii]
  end
  return x
end


text:
"""
This file is for creating RFF structure and for fast evaluations at given inputs.
"""

using Random

# A set of random fourier functions
struct RandomFourierFunctions{T<:AbstractFloat}
	directions::Matrix{T} # random direction
    offsets::Vector{T} # random shift in [0,2*pi]
    sigma::T

    function RandomFourierFunctions{T}(bandwidth, input, num_functions) where T<:AbstractFloat
        w = randn(T, num_functions, input) / T(bandwidth)
        # draw random offsets from a uniform distribution in [-pi,pi]
        b = rand(T, (num_functions,)) * T(2pi) .- T(pi)
        new{T}(w, b, T(bandwidth))
    end
end

function (rff::RandomFourierFunctions)(x::AbstractArray)
    W, b = rff.directions, rff.offsets
	nfeat = size(W,1)
    sD = eltype(b)(sqrt(2/nfeat))
    sD .* cos.(W*x .+ b)
end
(rff::RandomFourierFunctions{T})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}} = invoke(rff, Tuple{AbstractArray}, x)
(rff::RandomFourierFunctions{T})(x::AbstractArray{<:AbstractFloat}) where {T <: Union{Float32,Float64}} = rff(T.(x))

function (rff::RandomFourierFunctions{T})(z::AbstractArray{T, N},
                                          x::AbstractArray{T, N}) where {T <: Union{Float32,Float64}, N}
    nfeat = size(rff.directions,1)
    sD = T(sqrt(2/nfeat))

    W, b = rff.directions, rff.offsets
    mul!(z, W, x)
    z .+= b
    z .= sD .* cos.(z)
    z
end


text:
function __testset_nouncertainty_shared(r::RoutingSolution, rd::RoutingData, dm)
    @test r !== nothing
    @test r.data == rd
    @test r.n_iter == 1
    @test r.n_matrices == 1
    @test r.n_cuts == 0
    @test r.n_columns == 0
    @test r.time_precompute_ms > 0.0
    @test r.total_time_solve_ms > 0.0
    @test r.total_time_export_ms == 0.0
    @test r.total_time_ms > 0.0
    @test length(r.objectives) == 1
    @test length(r.matrices) == 1
    @test length(r.routings) == 1
    @test r.master_model !== nothing
    @test r.matrices[1][1] == dm
end

function __testset_nouncertainty_minmax(edge_obj, type, cg, opt, g, paths, k, d, dm)
    @testset "Objective aggregation: minimum maximum" begin
        mt = ModelType(edge_obj, MinimumMaximum(), type, cg, Automatic(), NoUncertaintyHandling(), NoUncertainty())
        rd = RoutingData(g, k, opt, mt, traffic_matrix=dm)

        r = nothing
        if typeof(edge_obj) == AlphaFairness && edge_obj.α != 0.0
            msg = "The objective function $(edge_obj) does not support minimisation. Proceed with caution."
            @test_logs (:warn, msg) r = compute_routing(rd)
        else
            @test_nowarn r = compute_routing(rd)
        end

        __testset_nouncertainty_shared(r, rd, dm)

        status = termination_status(r.master_model.model)
        if edge_obj == Load()
            @test status == MOI.OPTIMAL
            @test r.objectives[1] ≈ 0.571428 atol=1.0e-5
        elseif edge_obj == KleinrockLoad()
            @test status == MOI.OPTIMAL
            @test r.objectives[1] ≈ 1.333333 atol=1.0e-5
        elseif edge_obj == FortzThorupLoad()
            @test status == MOI.OPTIMAL
            @test r.objectives[1] ≈ 1.047619 atol=1.0e-5
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 0.0
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] ≈ 0.571428 atol=1.0e-5
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α in [0.5, 1.0, 1.5, 2.0]
            # Minimising the total fairness, which is usually nonsensical
            # (hence no objective value check).
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
        else
            println(r.objectives)
            @test_broken true == false
        end

        if status == MOI.OPTIMAL
            # Don't perform these tests for MOI.ALMOST_OPTIMAL: most constraint
            # are respected, but the sum over all paths may be "slightly" off
            # (1.05, for instance).

            @test r.routings[1].data == rd
            @test length(r.routings[1].demands) == ne(k) # Number of demands
            @test r.routings[1].demands[1] == d
            @test length(r.routings[1].path_flows) == ne(k) # Number of demands
            @test length(r.routings[1].path_flows[d]) in [2, 3] # Number of paths
            @test length(r.routings[1].edge_flows) == ne(k) # Number of demands
            @test length(r.routings[1].edge_flows[d]) in [3, 4, 5] # Number of edges

            total_flow = sum(x for x in values(r.routings[1].path_flows[d]))
            if abs(total_flow - 1.0) > 1.0e-3
                # SCS has troubles with some models (total_flow ≈ 0.998207).
                @test_broken total_flow ≈ 1.0 atol=1.0e-3
            else
                # Relaxed tolerance needed for some cases of α-fairness.
                @test total_flow ≈ 1.0 atol=1.0e-3
            end
        end
    end

    @testset "Objective aggregation: maximum minimum" begin
        mt = ModelType(edge_obj, MaximumMinimum(), type, cg, Automatic(), NoUncertaintyHandling(), NoUncertainty())
        rd = RoutingData(g, k, opt, mt, traffic_matrix=dm)

        r = nothing
        if typeof(edge_obj) in [KleinrockLoad, FortzThorupLoad]
            msg = "The objective function $(edge_obj) does not support maximisation. Proceed with caution."
            @test_logs (:warn, msg) r = compute_routing(rd)
        else
            @test_nowarn r = compute_routing(rd)
        end

        __testset_nouncertainty_shared(r, rd, dm)

        status = termination_status(r.master_model.model)
        if edge_obj == Load()
            @test status == MOI.OPTIMAL
            @test r.objectives[1] ≈ 0.571428 atol=1.0e-5
        elseif edge_obj == KleinrockLoad()
            @test status in [MOI.DUAL_INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED, MOI.SLOW_PROGRESS]
        elseif edge_obj == FortzThorupLoad()
            @test status in [MOI.DUAL_INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED, MOI.SLOW_PROGRESS]
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 0.0
            # Exactly as if optimising the load! Except that the power-cone model
            # may pose numerical problems.
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] ≈ 0.571428 atol=1.0e-5
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 0.5
            # SCS has troubles with this one: it sometimes thinks it is
            # unbounded, othertimes it gives a wrong answer.
            if status in [MOI.ALMOST_DUAL_INFEASIBLE, MOI.DUAL_INFEASIBLE]
                @test_broken status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
                @test_broken r.objectives[1] ≈ 1.511857 atol=1.0e-5
            else
                @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
                if status == MOI.OPTIMAL
                    if abs(r.objectives[1] - 1.511857) < 1.0e-5
                        @test r.objectives[1] ≈ 1.511857 atol=1.0e-5
                    else
                        @test_broken r.objectives[1] ≈ 1.511857 atol=1.0e-5
                    end
                end
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 1.0
            @test status == MOI.OPTIMAL
            @test r.objectives[1] ≈ -0.559636 atol=1.0e-4 # Looks like there are slight variations betweens flows and paths.
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 1.5
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL && abs(r.objectives[1] + 2.645751) < 1.0e-5
                # Neither SCS nor ECOS can find the optimum.
                # CPLEX (SOCP): -1.322875
                # Mosek (SOCP and POW): -2.645751
                @test r.objectives[1] ≈ -2.645751 atol=1.0e-5
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 2.0
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                if abs(r.objectives[1] + 1.749999) < 1.0e-5
                    @test r.objectives[1] ≈ -1.749999 atol=1.0e-5
                else
                    # Neither SCS nor ECOS can find the optimum.
                    @test_broken r.objectives[1] ≈ -1.749999 atol=1.0e-5
                end
            end
        else
            println(r.objectives)
            @test_broken true == false
        end

        if status == MOI.OPTIMAL
            # Don't perform these tests for MOI.ALMOST_OPTIMAL: most constraint
            # are respected, but the sum over all paths may be "slightly" off
            # (1.05, for instance).

            @test r.routings[1].data == rd
            @test length(r.routings[1].demands) == ne(k) # Number of demands
            @test r.routings[1].demands[1] == d
            @test length(r.routings[1].path_flows) == ne(k) # Number of demands
            @test length(r.routings[1].path_flows[d]) in [2, 3] # Number of paths
            @test length(r.routings[1].edge_flows) == ne(k) # Number of demands
            @test length(r.routings[1].edge_flows[d]) in [3, 4, 5] # Number of edges
            @test sum(x for x in values(r.routings[1].path_flows[d])) ≈ 1.0 atol=1.0e-3
        end
    end
end

function __testset_nouncertainty_mintot(edge_obj, type, cg, opt, g, paths, k, d, dm)
    @testset "Objective aggregation: minimum total" begin
        mt = ModelType(edge_obj, MinimumTotal(), type, cg, Automatic(), NoUncertaintyHandling(), NoUncertainty())
        rd = RoutingData(g, k, opt, mt, traffic_matrix=dm)

        r = nothing
        if typeof(edge_obj) == AlphaFairness && edge_obj.α != 0.0
            msg = "The objective function $(edge_obj) does not support minimisation. Proceed with caution."
            @test_logs (:warn, msg) r = compute_routing(rd)
        else
            @test_nowarn r = compute_routing(rd)
        end
        @test r !== nothing

        __testset_nouncertainty_shared(r, rd, dm)

        status = termination_status(r.master_model.model)
        if edge_obj == Load() || (typeof(edge_obj) == AlphaFairness && edge_obj.α == 0.0)
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] ≈ 2.0 atol=1.0e-5
            end
        elseif edge_obj == KleinrockLoad()
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] ≈ 5.952135 atol=1.0e-5
            end
        elseif edge_obj == FortzThorupLoad()
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] ≈ 4.666666 atol=1.0e-5
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 0.5
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] ≈ 0.000000 atol=1.0e-5
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α in [0.5, 1.0, 1.5, 2.0]
            # Minimising the total fairness, which is usually nonsensical.
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] <= -1000.0
            end
        else
            println(r.objectives)
            @test_broken true == false
        end

        if status == MOI.OPTIMAL
            # Don't perform these tests for MOI.ALMOST_OPTIMAL: most constraint
            # are respected, but the sum over all paths may be "slightly" off
            # (1.05, for instance).

            # This solution only uses two paths out of three, i.e. three edges out of five.
            @test r.routings[1].data == rd
            @test length(r.routings[1].demands) == ne(k) # Number of demands
            @test r.routings[1].demands[1] == d
            @test length(r.routings[1].paths) in [2, 3] # Number of paths
            @test length(r.routings[1].path_flows) == ne(k) # Number of demands
            @test length(r.routings[1].path_flows[d]) in [2, 3] # Number of paths
            @test length(r.routings[1].edge_flows) == ne(k) # Number of demands
            @test length(r.routings[1].edge_flows[d]) in [3, 4, 5] # Number of edges

            total_flow = sum(x for x in values(r.routings[1].path_flows[d]))
            if abs(total_flow - 1.0) > 1.0e-3
                # SCS has troubles with this one (total_flow ≈ 0.977865 for
                # 2-fair, or even 0.998207 or 1.004190 for others).
                @test_broken total_flow ≈ 1.0 atol=1.0e-3
            else
                @test total_flow ≈ 1.0 atol=1.0e-3
            end
        end
    end

    @testset "Objective aggregation: maximum total" begin
        mt = ModelType(edge_obj, MaximumTotal(), type, cg, Automatic(), NoUncertaintyHandling(), NoUncertainty())
        rd = RoutingData(g, k, opt, mt, traffic_matrix=dm)
        r = nothing

        if typeof(edge_obj) != AlphaFairness && typeof(edge_obj) != Load
            msg = "The objective function $(edge_obj) does not support maximisation. Proceed with caution."
            @test_logs (:warn, msg) r = compute_routing(rd)
        else
            @test_nowarn r = compute_routing(rd)
        end

        __testset_nouncertainty_shared(r, rd, dm)

        status = termination_status(r.master_model.model)
        if edge_obj == Load()
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] ≈ 4.0 atol=1.0e-5
            end
        elseif edge_obj == KleinrockLoad() || edge_obj == FortzThorupLoad()
            # Maximising these objective functions does not really make sense.
            @test status in [MOI.INFEASIBLE_OR_UNBOUNDED, MOI.DUAL_INFEASIBLE]
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 0.0
            @test status == MOI.OPTIMAL
            @test r.objectives[1] ≈ 4.0 atol=1.0e-5
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 0.5
            if ! edge_obj.force_power_cone || (edge_obj.force_power_cone && status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL])
                # SCS incorrectly determines this problem is unbounded with
                # the power-cone formulation, but Mosek sees no problem
                # and yields the same solution as for the SOCP formulation.
                @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
                if status == MOI.OPTIMAL
                    @test r.objectives[1] ≈ 8.326663 atol=1.0e-5
                end
            else
                @test_broken status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
                @test_broken r.objectives[1] ≈ 8.326663 atol=1.0e-5
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 1.0
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                @test r.objectives[1] ≈ -2.214 atol=1.0e-3
                # Mosek: -2.214330
                # ECOS: -2.214426.
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 1.5
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                # SCS has numerical problems here: it outputs
                # -1.597694850698785e-7, with a solution that is not really
                # feasible.
                if abs(r.objectives[1] + 12.696022) <= 1.0e-5
                    @test r.objectives[1] ≈ -12.696022 atol=1.0e-5
                else
                    @test_broken r.objectives[1] ≈ -12.696022 atol=1.0e-5
                end
            end
        elseif typeof(edge_obj) == AlphaFairness && edge_obj.α == 2.0
            @test status in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL]
            if status == MOI.OPTIMAL
                if abs(r.objectives[1] + 8.214101) < 1.0e-5
                    @test r.objectives[1] ≈ -8.214101 atol=1.0e-5
                else
                    # Neither SCS nor ECOS can find the optimum.
                    @test_broken r.objectives[1] ≈ -8.214101 atol=1.0e-5
                end
            end
        else
            println(r.objectives)
            @test_broken true == false
        end

        if status == MOI.OPTIMAL
            # Don't perform these tests for MOI.ALMOST_OPTIMAL: most constraint
            # are respected, but the sum over all paths may be "slightly" off
            # (1.05, for instance).

            # This solution only uses two paths out of three, i.e. three edges out of five.
            @test r.routings[1].data == rd
            @test length(r.routings[1].demands) == ne(k) # Number of demands
            @test r.routings[1].demands[1] == d
            @test length(r.routings[1].paths) in [2, 3] # Number of paths
            @test length(r.routings[1].path_flows) == ne(k) # Number of demands
            @test length(r.routings[1].path_flows[d]) in [2, 3] # Number of paths
            @test length(r.routings[1].edge_flows) == ne(k) # Number of demands
            @test length(r.routings[1].edge_flows[d]) in [3, 4, 5] # Number of edges
            @test sum(x for x in values(r.routings[1].path_flows[d])) ≈ 1.0 atol=1.0e-3
        end
    end
end

function __testset_nouncertainty_mmf(edge_obj, type, cg, opt, g, paths, k, d, dm)
    @testset "Objective aggregation: min-max fair" begin
        mt = ModelType(edge_obj, MinMaxFair(0.001), type, cg, Automatic(), NoUncertaintyHandling(), NoUncertainty())
        rd = RoutingData(g, k, opt, mt, traffic_matrix=dm, logmessage=(msg)->nothing)
        r = compute_routing(rd)

        # Similar to __testset_nouncertainty_shared, but many things differ for
        # MMF. Mostly, it's an iterative process.
        @test r !== nothing
        @test r.data == rd
        @test r.n_iter in [1, 2]
        @test r.n_matrices == 1
        @test r.n_cuts == 0
        @test r.n_columns == 0
        @test r.time_precompute_ms > 0.0
        @test r.total_time_solve_ms > 0.0
        @test r.total_time_export_ms == 0.0
        @test r.total_time_ms > 0.0
        @test length(r.matrices) == 1
        @test length(r.routings) in [1, 2]
        @test length(r.objectives) in [1, 2]
        @test r.master_model !== nothing
        @test r.matrices[1][1] == dm

        @test r.routings[1].data == rd
        @test length(r.routings[1].demands) == ne(k) # Number of demands
        @test r.routings[1].demands[1] == d
        @test length(r.routings[1].paths) <= length(paths) # Number of paths
        @test length(r.routings[1].paths) >= 1
        @test length(r.routings[1].path_flows) == ne(k) # Number of demands
        @test length(r.routings[1].path_flows[d]) in [1, 2, 3] # Number of paths
        @test length(r.routings[1].edge_flows) == ne(k) # Number of demands
        @test length(r.routings[1].edge_flows[d]) in [1, 3, 5] # Number of edges
        @test sum(x for x in values(r.routings[1].path_flows[d])) ≈ 1.0 atol=1.0e-5

        # Check that all paths are used, i.e. each of them sees at least 20% of the traffic.
        @test all(collect(values(r.routings[1].path_flows[d])) .>= 0.2)

        if typeof(edge_obj) != AlphaFairness # Probably due to poor accuracy.
            @test collect(values(r.routings[1].path_flows[d])) ≈ [0.2857, 0.2857, 0.4285] atol=1.0e-4
        end
    end

    @testset "Objective aggregation: max-min fair" begin
        mt = ModelType(edge_obj, MaxMinFair(0.001), type, cg, Automatic(), NoUncertaintyHandling(), NoUncertainty())
        rd = RoutingData(g, k, opt, mt, traffic_matrix=dm, logmessage=(msg)->nothing)
        r = compute_routing(rd)

        # Similar to __testset_nouncertainty_shared, but many things differ for
        # MMF. Mostly, it's an iterative process.
        @test r !== nothing
        @test r.data == rd
        @test r.n_iter in [1, 2]
        @test r.n_matrices == 1
        @test r.n_cuts == 0
        @test r.n_columns == 0
        @test r.time_precompute_ms > 0.0
        @test r.total_time_solve_ms > 0.0
        @test r.total_time_export_ms == 0.0
        @test r.total_time_ms > 0.0
        @test length(r.matrices) == 1
        @test length(r.routings) in [1, 2]
        @test length(r.objectives) in [1, 2]
        @test r.master_model !== nothing
        @test r.matrices[1][1] == dm

        @test r.routings[1].data == rd
        @test length(r.routings[1].demands) == ne(k) # Number of demands
        @test r.routings[1].demands[1] == d
        @test length(r.routings[1].paths) <= length(paths) # Number of paths
        @test length(r.routings[1].paths) >= 1
        @test length(r.routings[1].path_flows) == ne(k) # Number of demands
        @test length(r.routings[1].path_flows[d]) in [1, 2, 3] # Number of paths
        @test length(r.routings[1].edge_flows) == ne(k) # Number of demands
        @test length(r.routings[1].edge_flows[d]) in [1, 3, 5] # Number of edges
        # TODO: is something broken?
        # @test sum(x for x in values(r.routings[1].path_flows[d])) ≈ 1.0 atol=1.0e-5
        #
        # # Check that all paths are used, i.e. each of them sees at least 20% of the traffic.
        # @test all(collect(values(r.routings[1].path_flows[d])) .>= 0.2)
        #
        # if typeof(edge_obj) != AlphaFairness # Probably due to poor accuracy.
        #     @test collect(values(r.routings[1].path_flows[d])) ≈ [0.2857, 0.2857, 0.4285] atol=1.0e-4
        # end
    end
end

# Required cones: SOCP (Kleinrock, α-fairness), EXP (α-fairness), POW (α-fairness).
# ECOS only provides SOCP and EXP, use SCS for the rest (worse accuracy that ECOS).
opt_ecos = ECOS.Optimizer
opt_scs = SCS.Optimizer

g = MetaDiGraph()
add_vertex!(g, :name, "A")
add_vertex!(g, :name, "B")
add_vertex!(g, :name, "C")
add_vertex!(g, :name, "D")
add_edge!(g, 1, 2, :capacity, 2.0)
add_edge!(g, 1, 3, :capacity, 2.0)
add_edge!(g, 2, 4, :capacity, 2.0)
add_edge!(g, 3, 4, :capacity, 2.0)
add_edge!(g, 1, 4, :capacity, 3.0)
paths = [[Edge(1, 4)], [Edge(1, 2), Edge(2, 4)], [Edge(1, 3), Edge(3, 4)]]

k = SimpleDiGraph(4)
add_edge!(k, 1, 4)
d = Edge(1, 4)
dm = Dict{Edge{Int}, Float64}(d => 4.0)

@testset "Formulation type: $type$(ifelse(cg, " with column generation", ""))" for (type, cg) in [(FlowFormulation(), false), (PathFormulation(), false), (PathFormulation(), true)]
    @testset "Edge-wise objective: $edge_obj" for edge_obj in [
            Load(), KleinrockLoad(), FortzThorupLoad(),
            AlphaFairness(0.0),
            AlphaFairness(0.5, false), AlphaFairness(0.5, true),
            AlphaFairness(1.0),
            AlphaFairness(1.5, true), AlphaFairness(1.5, false),
            AlphaFairness(2.0, true), AlphaFairness(2.0, false)
        ]
        if cg && edge_obj != Load()
            # No column generation for the nonlinear case.
            # TODO: study something like http://liu.diva-portal.org/smash/get/diva2:841508/FULLTEXT01.pdf?
            # TODO: or like https://www.jstor.org/stable/2628397?seq=1?
            continue
        end

        @testset "No uncertainty" begin
            opt = opt_ecos
            if typeof(edge_obj) == AlphaFairness && edge_obj.force_power_cone
                opt = opt_scs
            end
            # opt = Mosek.Optimizer
            # opt = CPLEX.Optimizer

            if ! cg # TODO: not yet implemented.
                __testset_nouncertainty_minmax(edge_obj, type, cg, opt, g, paths, k, d, dm)
                __testset_nouncertainty_mintot(edge_obj, type, cg, opt, g, paths, k, d, dm)
            end

            if typeof(edge_obj) != AlphaFairness || (! edge_obj.force_power_cone && edge_obj.α != 0.0)
                # Power-cone models cannot be solved to optimality with ECOS or SCS,
                # and they also sometimes have problems with linear models.
                if ! (typeof(edge_obj) == AlphaFairness && edge_obj.α >= 1.0) # TODO: MMF doesn't work there: a bug somewhere?
                    __testset_nouncertainty_mmf(edge_obj, type, cg, opt, g, paths, k, d, dm)
                end
            end
        end
    end
end


text:
first = 1
second = 2.0
stringvariable = "Hello, World!"
character = 'F'
boolval = true

println(first)
println(second)
println(stringvariable)
println(character)
println(boolval)

max_stars_repo_name: ryanorsinger/JuliaPlayground

max_issues_repo_name: ryanorsinger/JuliaPlayground

max_forks_repo_name: ryanorsinger/JuliaPlayground


text:
@testset "8 Special functions" begin
    include("8.1 Error functions.jl")
    include("8.10 Formal derivatives.jl")
    include("8.2 Fresnel integral functions.jl")
    include("8.3 Exponential integral functions.jl")
    include("8.4 Trig integral functions.jl")
    include("8.5 Hyperbolic integral functions.jl")
    include("8.6 Gamma functions.jl")
    include("8.7 Zeta function.jl")
    include("8.8 Polylogarithm function.jl")
    include("8.9 Product logarithm function.jl")
end

max_stars_repo_path: test/8 Special functions/.jl

max_issues_repo_path: test/8 Special functions/.jl

max_forks_repo_path: test/8 Special functions/.jl


text:
using Measures
include("helper_functions.jl")

begin
    λ1 = 1
    λ2 = 1.5
    α = 1
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    parameters = (λ1, λ2, α, μ, ν)
    opponent_parameters = (λ2, λ1, α, μ, ν)
    p, q = equilibrium_probability(parameters...)
    qs = 0:0.001:1
end

p1_cost = cost.(pstar.(qs, parameters...), qs, parameters...)
p2_cost = cost.(qs, pstar.(qs, parameters...), opponent_parameters...)
p2_mask = p2_cost .>= 0
p1_mask = p1_cost .>= 0

begin
    cost_plot = plot(qs[p1_mask], p1_cost[p1_mask], label = "Player 1 cost")
    plot!(qs[p2_mask], p2_cost[p2_mask], label = "Player 2 cost")
    vline!([q], label = "Equilibrium", lc = :red, ls = :dash)
    xlabel!("Player 2 probability")
    ylabel!("Cost")
    plot!(size = (600, 250), ylims = (5, 20), xlims = (0.57, 0.68), margin = 2.5mm)
    plot!(legend = :bottomright)
end
# savefig("outputs/costs_vary_q_15.pdf")

begin
    total_cost_var = p1_cost + p2_cost
    p1gp2 = sign.(p1_cost[p1_mask] - p2_cost[p1_mask])
    tc_plot = plot(cost_plot, qs, total_cost_var, label = "Total cost")
    plot!(xlims = (0.57, 0.65), ylims = (0, 40))
    mincost, pos = findmin(total_cost_var[0.56 .< qs .< 0.7])
    vline!([qs[0.56 .< qs .< 0.7][pos]], label = "Lower total cost", ls = :dash, lc = :black)
    vline!([qs[p1_mask][p1gp2 .!= circshift(p1gp2, 1)]], label = "Minmax cost", ls = :dash)
    plot!(size = (600, 300), legend = :topright)
end
savefig("outputs/soc_opt_15.pdf")

calculate_assumptions(0.5, 0.5, parameters...)

# function total_cost(p, q, λ1, λ2, α, μ, ν)
#     p1_cost = cost(p, q, λ1, λ2, α, μ, ν)
#     p2_cost = cost(q, p, λ2, λ1, α, μ, ν)

#     return (p1_cost, p2_cost, (p1_cost > 0 ? p1_cost : 100000) + (p2_cost > 0 ? p2_cost : 100000))
# end

# begin
#     qn = q - 0.004
#     pn = pstar(qn, parameters...)
#     total_cost(p, q, parameters...) .- total_cost(pn, qn, parameters...)
#     # println("q: ", q)
#     # println("qn: ", qn)
#     # println("p: ", p)
#     # println("pn: ", pn)
# end

total_cost(p, q, parameters...)
total_cost(p, 0.5, parameters...)
total_cost(pn, qn, parameters...)

res = total_cost.(qs, pstar.(qs, parameters...), parameters...)
p1_c = getindex.(res, 1)
p2_c = getindex.(res, 2)
tcs = getindex.(res, 3)
plot(qs, p1_c)
plot!(qs, p2_c)
plot!(qs, tcs)
plot!(ylims = (-50, 50))
vline!([pn])
vline!([q])

plot(
    0:0.005:1,
0:0.005:1,
(p, q) -> total_cost(p, q, parameters...),
    st = :surface,
    zlims = (-50, 50),
    camera = (90, 90))
scatter!([p], [q], [total_cost(p, q, parameters...)], legend = nothing)

2
# q
# perturbed_q = 0.6
# current_cost = cost(p, q, parameters...)
# p1_newcost = cost(pstar(perturbed_q, parameters...), perturbed_q, parameters...)
# analytic_best_responded_cost(q, parameters...)
# p2_newcost = cost(perturbed_q, pstar(perturbed_q, parameters...), opponent_parameters...)

# r2perturb = pstar(perturbed_q, parameters...)
# r3p = pstar(r2perturb, opponent_parameters...)
# r4p = pstar(r3p, parameters...)

# ps = pstar.(qs, parameters...)

# p1_plot = plot(qs, ps, label = "p*", lims = (0, 1))
# vline!(p1_plot, [perturbed_q], ls = :dash)
# vline!(p1_plot, [r2perturb], ls = :dash)
# vline!(p1_plot, [r3p], ls = :dash)
# vline!(p1_plot, [r4p], ls = :dash)

# p1_plot = plot(qs, pstar.(qs, opponent_parameters...), label = "q*", lims = (0, 1))
# vline!([q], label = "Pareto probability", lc = :red, ls = :dash)
# xlabel!("Opponent's probability")
# ylabel!("Cost")

# p1_queue_length = (1 .- ps) .* λ1 .* private_cost.(ps, qs, parameters...)
# p2_queue_length = (1 .- qs) .* λ2 .* private_cost.(qs, ps, opponent_parameters...)
# public_queue_length = ((λ1 .* ps) + (λ2 .* qs)) .* public_cost.(ps, qs, parameters...)

# plot(qs, p1_queue_length, label = "Player 1 queue length")
# plot!(qs[p2_queue_length .> 0], p2_queue_length[p2_queue_length .> 0], label = "Player 2 queue length")
# plot!(qs, public_queue_length, label = "Public queue length")
# ylims!(0,30)

# q
# pstar(q, parameters...)


text:
module BiobakeryUtilsTests

using Random
using ReTest
using BiobakeryUtils
using BiobakeryUtils.Conda
using SparseArrays
using DelimitedFiles
using CSV

isdir(Conda.bin_dir(:BiobakeryUtils)) || BiobakeryUtils.install_deps()
ENV["PATH"] = ENV["PATH"] * ":" * Conda.bin_dir(:BiobakeryUtils)

@testset "CLI" begin
    @testset "Utilities" begin
        cmd = ["thing", "foo_bar"]
        cmd2 = copy(cmd)

        BiobakeryUtils.add_cli_kwargs!(cmd, Dict(:some_thing=> "foo", :bool=> true))
        @test all(cmd .== ["thing", "foo_bar", "--some_thing", "foo", "--bool"])
        BiobakeryUtils.add_cli_kwargs!(cmd2, Dict(:some_thing=> "foo", :bool=> true); optunderscores=false)
        @test all(cmd2 .== ["thing", "foo_bar", "--some-thing", "foo", "--bool"])
    end

    @testset "Metaphlan" begin
        @test BiobakeryUtils.check_for_install("metaphlan") |> isnothing
        @test BiobakeryUtils.check_for_install("merge_metaphlan_tables.py") |> isnothing

        @test metaphlan("", ""; help=true).exitcode == 0

        profiles = filter(f-> contains(f, "_profile.tsv"), readdir(joinpath(@__DIR__, "files/metaphlan"), join=true))
        @test metaphlan_merge(profiles, joinpath(@__DIR__, "files/metaphlan/merged_abundance_table.tsv")).exitcode == 0
    end

    @testset "Humann" begin
        @test BiobakeryUtils.check_for_install("humann") |> isnothing
        @test BiobakeryUtils.check_for_install("humann_rename_table") |> isnothing
        @test BiobakeryUtils.check_for_install("humann_renorm_table") |> isnothing
        @test BiobakeryUtils.check_for_install("humann_join_tables") |> isnothing
        @test BiobakeryUtils.check_for_install("humann") |> isnothing
        @test humann("", ""; help=true).exitcode == 0

    end
end

@testset "Metaphlan" begin
    profile_1 = metaphlan_profile(joinpath(@__DIR__, "files/metaphlan/SRS014464-Anterior_nares_profile.tsv");  sample="SRS014464")
    @test profile_1["k__Bacteria", "SRS014464"] == 100.0
    @test profile_1["o__Pseudomonadales", "SRS014464"] == 97.28734
    @test size(profile_1) == (13, 1)
    profile_2 = metaphlan_profile(joinpath(@__DIR__, "files/metaphlan/SRS014459-Stool_profile.tsv"), 3)
    @test size(profile_2) == (2, 1)
    @test profile_2["p__Firmicutes", "SRS014459-Stool_profile"] == 68.90167
    profile_3 = metaphlan_profile(joinpath(@__DIR__, "files/metaphlan/SRS014464-Anterior_nares_profile.tsv"), :phylum)
    @test size(profile_3) == (2, 1)
    @test profile_3["p__Proteobacteria", 1] == 97.28734

    merge_profile_1 = metaphlan_profiles(joinpath(@__DIR__, "files/metaphlan/merged_abundance_table.tsv"); samplestart=3)
    @test size(merge_profile_1) == (62, 6)
    @test merge_profile_1["g__Moraxella", 5] == 97.28734
    merge_profile_2 = metaphlan_profiles(joinpath(@__DIR__, "files/metaphlan/merged_abundance_table.tsv"), :family; samplestart=3)
    @test size(merge_profile_2) == (13, 6)
    @test merge_profile_2["f__Micrococcaceae", "SRS014464-Anterior_nares"] == 0.0
    merge_profile_3 = metaphlan_profiles(joinpath(@__DIR__, "files/metaphlan/merged_abundance_table.tsv"), 7; samplestart=3)
    @test size(merge_profile_3) == (16, 6)
    @test merge_profile_3["s__Haemophilus_haemolyticus", 3] == 1.35528
    CSV.write(joinpath(@__DIR__, "files/metaphlan/merged_abundance_table2.csv"), merge_profile_1)

    profiles = filter(f-> contains(f, "_profile.tsv"), readdir(joinpath(@__DIR__, "files/metaphlan"), join=true))
    @test_throws ArgumentError metaphlan_profiles(profiles; samples = ["sample1"])
    multi_profile_1 = metaphlan_profiles(profiles; samples=["sample$i" for i in 1:length(profiles)])
    @test abundances(multi_profile_1) == abundances(metaphlan_profiles(profiles))
    @test size(multi_profile_1) == (62, 6)
    @test multi_profile_1["p__Firmicutes", "sample1"] == 68.90167
    multi_profile_2 = metaphlan_profiles(profiles, 3; samples=["sample$i" for i in 1:length(profiles)])
    @test abundances(multi_profile_2) == abundances(metaphlan_profiles(profiles, :class))
    @test size(multi_profile_2) == (6,6)
    @test multi_profile_2["c__Bacteroidia", "sample1"] == 31.09833

    taxstring = "k__Archaea|p__Euryarchaeota|c__Methanobacteria|o__Methanobacteriales|f__Methanobacteriaceae|g__Methanobrevibacter|s__Methanobrevibacter_smithii"
    taxa = parsetaxa(taxstring)
    @test length(taxa) == 7
    @test parsetaxon(taxstring, 1) == Taxon("Archaea", :kingdom)
    @test parsetaxon(taxstring, :family) == Taxon("Methanobacteriaceae", :family)
    @test parsetaxon(taxstring) == Taxon("Methanobrevibacter_smithii", :species)
    @test_throws ArgumentError parsetaxon(taxstring, 8)

    @test parsetaxon("k__Archaea|p__Euryarchaeota|c__Methanobacteria", 2) == Taxon("Euryarchaeota", :phylum)
    @test parsetaxon("k__Archaea|p__Euryarchaeota|c__Methanobacteria") == Taxon("Methanobacteria", :class)
end

@testset "HUMAnN" begin
    p1 = humann_profile(joinpath(@__DIR__, "files/humann/single_1.tsv"))
    p2 = humann_profile(joinpath(@__DIR__, "files/humann/single_2.tsv"))
    @test p1 isa CommunityProfile
    @test size(p1) == (560, 1)
    @test samplenames(p1) == ["single_1"]
    @test samplenames(humann_profile(joinpath(@__DIR__, "files/humann/single_1.tsv"); sample = "sample1")) == ["sample1"]
    @test samplenames(humann_profile(joinpath(@__DIR__, "files/humann/single_1.tsv"); sample = MicrobiomeSample("sample1"))) == ["sample1"]

    @test all(f-> !hastaxon(f), features(p1)) # unstratified
    @test all(f-> !occursin('|', name(f)), features(p1))

    pj = humann_profiles(joinpath(@__DIR__, "files/humann/joined.tsv"))
    @test size(pj) == (560, 2)
    @test isempty(setdiff(features(pj), features(commjoin(p1, p2))))
    @test samplenames(pj) == samplenames(commjoin(p1, p2))

    pj_strat = humann_profiles(joinpath(@__DIR__, "files/humann/joined.tsv"); stratified = true)
    @test size(pj_strat) == (1358, 2)
    @test !isempty(setdiff(features(pj_strat), features(pj)))
    @test isempty(setdiff(featurenames(pj_strat), featurenames(pj)))
    @test isempty(setdiff(features(filter(!hastaxon, pj_strat)), features(pj)))
    CSV.write(joinpath(@__DIR__, "files/humann/joined_roundtrip.tsv"), pj_strat; delim='\t')
end

end # module

max_stars_repo_name: BioJulia/BiobakeryUtils

max_issues_repo_name: BioJulia/BiobakeryUtils

max_forks_repo_name: BioJulia/BiobakeryUtils


text:
using StaticArrays
function ZigZagBoomerang.move_forward!(τ, t, x::SVector, θ, Z::Union{BouncyParticle, ZigZag})
    t += τ
    x += θ .* τ
    t, x, θ
end
ZigZagBoomerang.smove_forward!(τ, t, x::SVector, θ, f, Z::Union{BouncyParticle, ZigZag}) = ZigZagBoomerang.move_forward!(τ, t, x, θ, Z)
function ZigZagBoomerang.reflect!(∇ϕx, x::SVector, θ, ::Boomerang)
    θ -= (2*dot(∇ϕx, θ)/ZigZagBoomerang.normsq(∇ϕx))*∇ϕx
    θ
end
function ZigZagBoomerang.reflect!(∇ϕx, x::SVector, θ, F::BouncyParticle)
    θ -= (2*dot(∇ϕx, θ)/ZigZagBoomerang.normsq(F.L\∇ϕx))*(F.L'\(F.L\∇ϕx))
    θ
end
function ZigZagBoomerang.refresh!(rng, θ::SVector{d}, F) where {d}
    ρ̄ = sqrt(1-F.ρ^2)
    F.ρ*θ + ρ̄*randn(rng, SVector{d})
end

function ZigZagBoomerang.refresh!(rng, θ::SVector{d}, F::BouncyParticle) where {d}
    ρ̄ = sqrt(1-F.ρ^2)
    F.ρ*θ + ρ̄*convert(typeof(θ), F.L'\randn(rng, SVector{d}))
end


text:
function plot_lindo(
    U;
    label="",
    limits = FRect(0, -(size(U,3)+1), size(U,2)+1, (size(U,3)+1)),
    kwargs...)

    scatter(reshape(U[3,:,:], size(U,2)*size(U,3)), reshape(U[4,:,:], size(U,2)*size(U,3)); label=label, limits=limits, kwargs...)
end

function render(
    sol;
    fps=24,
    path = "$(pwd())/img/$(size(sol.u[1], 2))x$(size(sol.u[1], 3)).mp4",
    limits = FRect(0, -(size(sol.u[1],3)+1), size(sol.u[1],2)+1, (size(sol.u[1],3)+1)))

    tmax = sol.t[end]

    scene = plot_lindo(sol(0.0), limits = limits)
    record(scene, path, framerate=fps) do io
        for t = 0.0:1.0/fps:tmax
            u = sol(t)
            push!(scene.plots[end].input_args[1], reshape(u[3,:,:], size(u,2)*size(u,3)))
            push!(scene.plots[end].input_args[2], reshape(u[4,:,:], size(u,2)*size(u,3)))
            recordframe!(io)
            Base.CoreLogging.@logmsg(-1,"Rendering...",_id = :Render, message="$(t)s", progress=t/tmax)
        end
    end
    Base.CoreLogging.@logmsg(-1,"Rendering...",_id = :Render, progress="done")
    @info "Saved animation to $path"
end


text:
# Visualization
function visualize!(mvis, model::KukaParticle, q;
		verbose = false, r_ball = 0.025, Δt = 0.1)

	setobject!(vis["ball"], Sphere(Point3f0(0),
				convert(Float32, r_ball)),
				MeshPhongMaterial(color = RGBA(51.0 / 255.0, 1.0, 1.0, 1.0)))

	settransform!(vis["ball"], compose(Translation(0.66,3.0,0.0)))

	state = model.state_cache1[Float64]

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	T = length(q)
    for t = 1:T
        q_kuka = kuka_q(q[t])
		q_particle = particle_q(q[t])
		set_configuration!(state, kuka_q(q[t]))

        MeshCat.atframe(anim,t) do
			set_configuration!(mvis,q_kuka)
            settransform!(vis["ball"], compose(Translation(q_particle), LinearMap(RotZ(0))))
		end
    end

    MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Primes

println(sum(primes(2000000)))

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
## Exercise 8 - 7
## Read the documentation of the string functions at https://docs.julialang.org/en/v1/manual/strings/. You might want to experiment with some of them to make sure you understand how they work. strip and replace are particularly useful.
## The documentation uses a syntax that might be confusing. For example, in search(string::AbstractString, chars::Chars, [start::Integer]), the brackets indicate optional arguments. So string and chars are required, but start is optional.
println("Ans: ")

myvar = "this is from myvar"
mynum = 10
println("'$myvar' and my num is $mynum")

println("Please enter something: ");
mystr = strip(readline())
println(replace(mystr, "g" => "f"))

println(firstindex(mystr, "g"))
println("End.")


text:
using Empirikos
using Test

discr = DiscretePriorClass(Empirikos.DataBasedDefault())

Zs = StandardNormalSample.([1.0;2.0])
discr_instantiated = Empirikos.set_defaults(discr, Zs)
@test extrema(discr_instantiated.support) == (1.0 - 1e-4, 2.0 + 1e-4)
@test length(discr_instantiated.support) == 300

discr_instantiated_hinted = Empirikos.set_defaults(discr, Zs; hints=:prior_grid_length=>991)
@test length(discr_instantiated_hinted.support) == 991
@test extrema(discr_instantiated_hinted.support) == extrema(discr_instantiated.support)

tmp_grid = 0:0.1:10
discr_grid = DiscretePriorClass(tmp_grid)
@test Empirikos.set_defaults(discr_grid, Zs) == discr_grid

npmle = NPMLE(discr, nothing)

npmle_instantiated = Empirikos.set_defaults(npmle, Zs)
@test npmle_instantiated.convexclass == discr_instantiated

npmle_typefield = NPMLE(discr, Float64)
npmle_typefield_instantiated = Empirikos.set_defaults(npmle_typefield, Zs)
@test npmle_instantiated.convexclass == npmle_typefield_instantiated.convexclass

npmle_grid = NPMLE(discr_grid, nothing)
@test Empirikos.set_defaults(npmle_grid, Zs) == npmle_grid


text:
using Knet
import Gym
import AutoGrad: value
import Random
using Statistics

# using TimerOutputs
# reset_timer!()

const F = Float32

mutable struct History
    xsize
    nA::Int
    γ::F
    states  # flattens all states into a vector
    actions::Vector{Int}
    rewards::Vector{F}
end

History(xsize, nA, γ, atype) = History(xsize, nA, γ,
                                convert(atype, zeros(F,0)),
                                zeros(Int, 0),
                                zeros(F,0))

function Base.push!(history, s, a, r)
    history.states = [history.states; s]
    push!(history.actions, a)
    push!(history.rewards, r)
end

function preprocess(I, atype)
    I = I[36:195,:,:]
    I = I[1:2:end, 1:2:end, 1]
    I[I .== 144] .= 0
    I[I .== 109] .= 0
    I[I .!= 0] .= 1
    return convert(atype, vec(I))
end

function discount(rewards, γ)
    R = similar(rewards)
    R[end] = rewards[end]
    for k = length(rewards)-1:-1:1
        R[k] = γ * R[k+1] + rewards[k]
    end
    # return R
    return (R .- mean(R)) ./ (std(R) + F(1e-10)) #speeds up training a lot
end

function initweights(atype, xsize, ysize)
    @assert xsize == (80, 80, 1)
    w = []
    ch1 = 4
    ch2 = 8
    nh = 64
    push!(w, xavier(8, 8, xsize[3], ch1))
    push!(w, zeros(1, 1, ch1, 1))
    push!(w, xavier(4, 4, ch1, ch2))
    push!(w, zeros(1, 1, ch2, 1))

    push!(w, xavier(nh, 9*9*ch2))
    push!(w, zeros(nh, 1))

    push!(w, xavier(ysize, nh)) # policy
    push!(w, zeros(ysize, 1))
    push!(w, xavier(1, nh)) # value
    push!(w, zeros(1, 1))

    return map(wi->convert(atype,wi), w)
end

function predict(w, x)
    x = reshape(x, 80, 80, 1, :)
    x = relu.(conv4(w[1], x, stride=4, padding=2) .+ w[2])
    x = relu.(conv4(w[3], x, stride=2) .+ w[4])
    x = mat(x)
    x = relu.(w[5]*x .+ w[6])
    prob_act = w[7] * x .+ w[8]
    value = w[9] * x .+ w[10]
    return prob_act, value
end

function sample_action(probs)
    probs = Array(probs)
    cprobs = cumsum(probs, dims=1)
    sampled = cprobs .> rand()
    actions = mapslices(argmax, sampled, dims=1)
    return actions[1]
end

function loss(w, history)
    xsize, nA = history.xsize, history.nA
    nS = prod(xsize)
    M = length(history.states) ÷ nS
    states = reshape(history.states, nS, M)

    p, V = predict(w, states)
    V = vec(V)

    R = discount(history.rewards, history.γ)
    R = convert(typeof(value(V)), R)
    A = R .- V   # advantage

    inds = history.actions + nA*(0:M-1)
    lp = logsoftmax(p, dims=1)[inds] # lp is a vector

    return -mean(lp .* value(A)) + L2Reg(A)
end

L2Reg(x) = mean(x .* x)

function main(;
        lr = 1e-2, # learning rate
        γ = 0.99, # discount rate
        episodes = 10000,  # max episodes played
        render = false, # if true display an episode every infotime ones
        seed = -1,
        infotime = 50,
        atype = gpu() >= 0 ? KnetArray{F} : Array{F},
    )

    @info "using $atype"
    env = Gym.GymEnv("Pong-v0")
    seed > 0 && (Random.seed!(seed); Gym.seed!(env, seed))
    xsize = (80, 80, 1) # input dimensions of the policy network
    nA = 3 # number of actions. Could also set to 2
    amap = Dict(1=>2, 2=>3, 3=>0) #up, down, no action

    w = initweights(atype, xsize, nA)
    opt = [Adam(lr=lr) for _=1:length(w)]
    avgreward = -21.0
    for episode=1:episodes
        state = Gym.reset!(env)
        episode_reward = 0
        history = History(xsize, nA, γ, atype)
        prev_s = convert(atype, zeros(prod(xsize)))
        while true
            curr_s = preprocess(state, atype)
            s = curr_s .- prev_s
            prev_s = curr_s

            p, V = predict(w, s)
            p = softmax(p, dims=1)
            action = sample_action(p)

            state, reward, done, info = Gym.step!(env, amap[action])

            push!(history, s, action, reward)
            episode_reward += reward
            episode % infotime == 0 && render && Gym.render(env)
            done && break
        end

        avgreward = 0.1*episode_reward + 0.9*avgreward
        println("Episode:$episode, Reward:$episode_reward, AvgReward:$(avgreward)")
        episode % infotime == 0 && Gym.close!(env)

        dw = grad(loss)(w, history)
        update!(w, dw, opt)
    end

    return w
end


text:
using KernelAbstractions, Test, CUDAapi
if CUDAapi.has_cuda_gpu()
    using CuArrays
    CuArrays.allowscalar(false)
end

# Simple kernel for matrix multiplication
@kernel function matmul_kernel!(a, b, c)
    i, j = @index(Global, NTuple)

    # creating a temporary sum variable for matrix multiplication
    tmp_sum = zero(eltype(c))
    for k = 1:size(a)[2]
        tmp_sum += a[i,k] * b[k, j]
    end

    c[i,j] = tmp_sum
end

# Creating a wrapper kernel for launching with error checks
function matmul!(a, b, c)
    if size(a)[2] != size(b)[1]
        println("Matrix size mismatch!")
        return nothing
    end
    if isa(a, Array)
        kernel! = matmul_kernel!(CPU(),4)
    else
        kernel! = matmul_kernel!(CUDA(),256)
    end
    kernel!(a, b, c, ndrange=size(c))
end

a = rand(256,123)
b = rand(123, 45)
c = zeros(256, 45)

# beginning CPU tests, returns event
ev = matmul!(a,b,c)
wait(ev)

@test isapprox(c, a*b)

# beginning GPU tests
if has_cuda_gpu()
    d_a = CuArray(a)
    d_b = CuArray(b)
    d_c = CuArray(c)

    ev = matmul!(d_a, d_b, d_c)
    wait(ev)

    @test isapprox(Array(d_c), a*b)
end


text:
export @settable
using MacroTools: prewalk
using MacroTools: splitdef, combinedef
using MacroTools: splitstructdef, isstructdef, combinestructdef

macro settable(ex)
    esc(settable(__module__, ex))
end

function arg_type(ex)::Tuple
    @match ex begin
        (name_::T_ = default_) => (name, T  )
        (name_ = default_    ) => (name, Any)
        (name_::T_           ) => (name, T  )
        (name_               ) => (name, Any)
    end
end

argsymbol(arg)::Symbol = first(arg_type(arg))
function argsymbol_typed(arg)::Expr
    name, T, = arg_type(arg)
    MacroTools.combinearg(name,T,false,nothing)
end

function has_posonly_constructor(dtype::Dict)
    fields = map(first, dtype[:fields])
    for constructor in dtype[:constructors]
        args = map(argsymbol, splitdef(constructor)[:args])
        if args == fields
            return true
        end
    end
    return false
end

function constructor_has_desired_fields(fields, dconstr::Dict)
    args = dconstr[:args]
    nargs = length(args)
    map(argsymbol, args) == fields[1:nargs] || return false
    Set(map(argsymbol, dconstr[:kwargs])) == Set(fields[nargs+1:end])
end

function best_constructor_template(dtype::Dict)::Dict
    fields = map(first, dtype[:fields])
    constructors = map(splitdef, dtype[:constructors])
    constructors = filter(d -> constructor_has_desired_fields(fields, d), constructors)
    if isempty(constructors)
        error("""
              There is no appropriate inner constructor.  At least one
              constructor has to have positional or keyword arguments with
              name matching with the struct fields.
              """)
    end
    nparams(c) = haskey(c, :params) ? length(c[:params]) : 0
    nargs(c) = haskey(c, :args) ? length(c[:args]) : 0
    min_nparams = minimum(nparams,constructors)
    constructors = filter(c -> nparams(c) == min_nparams, constructors)
    max_nargs = maximum(nargs, constructors)
    filter(c->nargs(c) == max_nargs, constructors)
    first(constructors)
end

function posonly_constructor_dict(dtype::Dict)::Dict
    def = best_constructor_template(dtype)::Dict
    fields = map(first, dtype[:fields])
    arg_dict = Dict(argsymbol(arg) => arg for arg in def[:args])
    kwarg_dict = Dict(argsymbol(arg) => arg for arg in def[:kwargs])
    newargs = []
    for field in fields
        if haskey(arg_dict, field)
            newarg = arg_dict[field]
        else
            newarg = kwarg_dict[field]
            delete!(kwarg_dict, field)
        end
        push!(newargs, argsymbol_typed(newarg))
    end
    newkwargs = collect(values(kwarg_dict))
    def[:args] = newargs
    def[:kwargs] = newkwargs
    def
end

function posonly_constructor(dtype::Dict)::Expr
    combinedef(posonly_constructor_dict(dtype))
end

function add_posonly_constructor(ex::Expr)::Expr
    dtype = splitstructdef(ex)
    if isempty(dtype[:constructors])
        ex
    elseif has_posonly_constructor(dtype)
        ex
    else
        push!(dtype[:constructors], posonly_constructor(dtype))
        @assert has_posonly_constructor(dtype)
        combinestructdef(dtype)
    end
end

function settable(M, code)::Expr
    code = macroexpand(M, code)
    MacroTools.postwalk(code) do ex
        ret = if isstructdef(ex)
            add_posonly_constructor(ex)
        else
            ex
        end
    end
end


text:
using MacroTools

function Base.convert(::Type{Expression}, ex::Expr)
    ex.head === :call || throw(ArgumentError("internal representation does not support non-call Expr"))

    op = eval(ex.args[1])  # HACK
    args = convert.(Expression, ex.args[2:end])

    return Operation(op, args)
end
Base.convert(::Type{Expression}, x::Expression) = x
Base.convert(::Type{Expression}, x::Number) = Constant(x)

function build_expr(head::Symbol, args)
    ex = Expr(head)
    append!(ex.args, args)
    ex
end

# used in parsing
isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
function flatten_expr!(x)
    isb = isblock(x)
    if isb
        x = MacroTools.striplines(x[1])
        filter!(z->z isa Symbol || z.head != :line, x.args)
        x = (x.args...,)
    end
    x
end

function build_function(rhss, vs, ps, args = (); version::FunctionVersion)
    var_pairs   = [(u.name, :(u[$i])) for (i, u) ∈ enumerate(vs)]
    param_pairs = [(p.name, :(p[$i])) for (i, p) ∈ enumerate(ps)]
    (ls, rs) = zip(var_pairs..., param_pairs...)

    var_eqs = Expr(:(=), build_expr(:tuple, ls), build_expr(:tuple, rs))

    if version === ArrayFunction
        X = gensym()
        sys_exprs = [:($X[$i] = $(convert(Expr, rhs))) for (i, rhs) ∈ enumerate(rhss)]
        let_expr = Expr(:let, var_eqs, build_expr(:block, sys_exprs))
        :(($X,u,p,$(args...)) -> $let_expr)
    elseif version === SArrayFunction
        sys_expr = build_expr(:tuple, [convert(Expr, rhs) for rhs ∈ rhss])
        let_expr = Expr(:let, var_eqs, sys_expr)
        :((u,p,$(args...)) -> begin
            X = $let_expr
            T = StaticArrays.similar_type(typeof(u), eltype(X))
            T(X)
        end)
    end
end

is_constant(::Constant) = true
is_constant(::Any) = false

is_operation(::Operation) = true
is_operation(::Any) = false

is_derivative(O::Operation) = isa(O.op, Differential)
is_derivative(::Any) = false

has_dependent(t::Variable) = Base.Fix2(has_dependent, t)
has_dependent(x::Variable, t::Variable) =
    any(isequal(t), x.dependents) || any(has_dependent(t), x.dependents)


text:
n = rand(1:38)
mat = matrixdepot(n)
matlist = matrixdepot(1:n)


text:
# doesn't fully work WIP
macro estrin(t, c...)
    N = length(c) ÷ 2

    tl = Expr(:tuple)
    tr = Expr(:tuple)
    for i = 1:N
        push!(tl.args, :(p[2($i)]))
        push!(tr.args, :(p[2($i)-1]))
    end

    a1 = gensym("a1")
    a2 = gensym("a2")
    a = gensym("a")

    blk = quote
        x = $(esc(t))
        p = $(Expr(:tuple, map(esc, c)...))

        $(:(@inbounds $(a1) = Vec{$N,Float64}($tl)))
        $(:(@inbounds $(a2) = Vec{$N,Float64}($tr)))

        # $(a1) = Vec{$N,Float64}($tl)
        # $(a2) = Vec{$N,Float64}($tr)
        $(a) =  muladd($(a1), x, $(a2))
    end

   while N > 1
        ext = Expr(:block, :(x = x * x))
        N = N÷2

        tl = Expr(:tuple)
        tr = Expr(:tuple)
        for i = 1:N
            push!(tl.args, :($a[2($i)]))
            push!(tr.args, :($a[2($i)-1]))
        end

        a1 = gensym("a1")
        a2 = gensym("a2")
        a = gensym("a")

        push!(ext.args, :(@inbounds $(a1) = Vec{$N,Float64}($tl)))
        push!(ext.args, :(@inbounds $(a2) = Vec{$N,Float64}($tr)))
        push!(ext.args, :($(a) = muladd($(a1), x, $(a2)) ))

        push!(blk.args,ext)
    end
    v = gensym("x")
    push!(blk.args,:(@inbounds $v =  $(a)[1]))
    push!(blk.args,:($v))
    return blk
end


text:
# ---
# title: 1342. Number of Steps to Reduce a Number to Zero
# id: problem1342
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/number-of-steps-to-reduce-a-number-to-zero/description/>
# hidden: true
# ---
#
# Given a non-negative integer `num`, return the number of steps to reduce it to
# zero. If the current number is even, you have to divide it by 2, otherwise,
# you have to subtract 1 from it.
#
#
#
# **Example 1:**
#
#
#
#     Input: num = 14
#     Output: 6
#     Explanation:
#     Step 1) 14 is even; divide by 2 and obtain 7.
#     Step 2) 7 is odd; subtract 1 and obtain 6.
#     Step 3) 6 is even; divide by 2 and obtain 3.
#     Step 4) 3 is odd; subtract 1 and obtain 2.
#     Step 5) 2 is even; divide by 2 and obtain 1.
#     Step 6) 1 is odd; subtract 1 and obtain 0.
#
#
# **Example 2:**
#
#
#
#     Input: num = 8
#     Output: 4
#     Explanation:
#     Step 1) 8 is even; divide by 2 and obtain 4.
#     Step 2) 4 is even; divide by 2 and obtain 2.
#     Step 3) 2 is even; divide by 2 and obtain 1.
#     Step 4) 1 is odd; subtract 1 and obtain 0.
#
#
# **Example 3:**
#
#
#
#     Input: num = 123
#     Output: 12
#
#
#
#
# **Constraints:**
#
#   * `0 <= num <= 10^6`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end