text:
# Optimization and solution query
#
#
# OBS: in xpress all functions with LP in name
#      actually refer to CONTINUOUS optimization
#      LP, QP, QCQP, SOCP
#      as opposed to Integer programming

"""
    optimize(model::Model)

Solve model
"""
function optimize(model::Model)
    @assert model.ptr_model != C_NULL
    if is_mip(model)
        mipoptimize(model)
    else
        lpoptimize(model)
    end
    nothing
end

"""
    lpoptimize(model::Model)

Solve models ignoring integrality.
Quadratic terms are NOT ignored.
"""
function lpoptimize(model::Model, flags::String="")
    @assert model.ptr_model != C_NULL
    start_time = time()
    ret = @xprs_ccall(lpoptimize, Cint, (Ptr{Nothing},Ptr{UInt8}),
        model.ptr_model, flags)
    model.time = time() - start_time
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

function mipoptimize(model::Model, flags::String="")
    @assert model.ptr_model != C_NULL
    start_time = time()
    ret = @xprs_ccall(mipoptimize, Cint, (Ptr{Nothing},Ptr{UInt8}),
        model.ptr_model, flags)
        model.time = time() - start_time
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

"""
    computeIIS(model::Model)

compute all ISS
"""
function computeIIS(model::Model)
    @assert model.ptr_model != C_NULL
    ret = @xprs_ccall(iisall, Cint, (Ptr{Nothing},), model.ptr_model)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

function loaddelayedrows(model::Model, mrows::Array{I}) where I<:Integer
    nrows = length(mrows)
    @assert model.ptr_model != C_NULL
    ret = @xprs_ccall(loaddelayedrows, Cint, (Ptr{Nothing},Cint,Ptr{Cint}),
        model.ptr_model,nrows,convert(Vector{Cint},mrows-1))
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

#################################################
#
#  solution status and optimization info
#
#################################################

# Unstarted (XPRS_LP_UNSTARTED).
# Optimal (XPRS_LP_OPTIMAL).
# Infeasible (XPRS_LP_INFEAS).
# Objective worse than cutoff (XPRS_LP_CUTOFF).
# Unfinished (XPRS_LP_UNFINISHED).
# Unbounded (XPRS_LP_UNBOUNDED).
# Cutoff in dual (XPRS_LP_CUTOFF_IN_DUAL).
# Problem is unsolved (XPRS_LP_UNSOLVED).
# Problem contains quadratic data, which is not convex (XPRS_LP_NONCONVEX).

@enum LP_Status LP_Unstarted LP_Optimal LP_Infeasible LP_CutOff LP_Unfinished LP_Unbounded LP_CutOffInDual LP_Unsolved LP_NonConvex LP_Unknown

function get_lp_status2(model::Model)
    code = get_lp_status_code(model)
    if code == XPRS_LP_UNSTARTED
        return LP_Unstarted
    elseif code == XPRS_LP_OPTIMAL
        return LP_Optimal
    elseif code == XPRS_LP_INFEAS
        return LP_Infeasible
    elseif code == XPRS_LP_CUTOFF
        return LP_CutOff
    elseif code == XPRS_LP_UNFINISHED
        return LP_Unfinished
    elseif code == XPRS_LP_UNBOUNDED
        return LP_Unbounded
    elseif code == XPRS_LP_CUTOFF_IN_DUAL
        return LP_CutOffInDual
    elseif code == XPRS_LP_UNSOLVED
        return LP_Unsolved
    elseif code == XPRS_LP_NONCONVEX
        return LP_NonConvex
    end
    return LP_Unknown
end

# XPRS_MIP_NOT_LOADED Problem has not been loaded.
# XPRS_MIP_LP_NOT_OPTIMAL Global search incomplete - the initial continuous relaxation has not been solved and no integer solution has been found.
# XPRS_MIP_LP_OPTIMAL Global search incomplete - the initial continuous relaxation has been solved and no integer solution has been found.
# XPRS_MIP_NO_SOL_FOUND Global search incomplete - no integer solution found.
# XPRS_MIP_SOLUTION Global search incomplete - an integer solution has been found.
# XPRS_MIP_INFEAS Global search complete - no integer solution found.
# XPRS_MIP_OPTIMAL Global search complete - integer solution found.
# XPRS_MIP_UNBOUNDED Global search incomplete - the initial continuous relaxation was found to be unbounded. A solution may have been found.

@enum MIP_Status MIP_NotLoaded MIP_LPNotOptimal MIP_LPOptimal MIP_NoSolFound MIP_Solution MIP_Infeasible MIP_Optimal MIP_Unbounded MIP_Unknown

function get_mip_status2(model::Model)
    code = get_mip_status_code(model)
    if code == XPRS_MIP_NOT_LOADED
        return MIP_NotLoaded
    elseif code == XPRS_MIP_LP_NOT_OPTIMAL
        return MIP_LPNotOptimal
    elseif code == XPRS_MIP_LP_OPTIMAL
        return MIP_LPOptimal
    elseif code == XPRS_MIP_NO_SOL_FOUND
        return MIP_NoSolFound
    elseif code == XPRS_MIP_SOLUTION
        return MIP_Solution
    elseif code == XPRS_MIP_INFEAS
        return MIP_Infeasible
    elseif code == XPRS_MIP_OPTIMAL
        return MIP_Optimal
    elseif code == XPRS_MIP_UNBOUNDED
        return MIP_Unbounded
    end
    return MIP_Unknown
end

const status_symbols_lp = Dict(
    XPRS_LP_UNSTARTED         => :unstarted,
    XPRS_LP_OPTIMAL           => :optimal,
    XPRS_LP_INFEAS            => :infeasible,
    XPRS_LP_CUTOFF            => :cutoff,
    XPRS_LP_UNFINISHED        => :unfinished,
    XPRS_LP_UNBOUNDED         => :unbounded,
    XPRS_LP_CUTOFF_IN_DUAL    => :cutoff_in_dual,
    XPRS_LP_UNSOLVED          => :unsolved,
    XPRS_LP_NONCONVEX         => :nonconvex
)

const status_symbols_mip = Dict(
    XPRS_MIP_NOT_LOADED       => :mip_not_loaded,
    XPRS_MIP_LP_NOT_OPTIMAL   => :mip_lp_not_optimal,
    XPRS_MIP_LP_OPTIMAL       => :mip_lp_optimal,
    XPRS_MIP_NO_SOL_FOUND     => :mip_no_sol_found,
    XPRS_MIP_SOLUTION         => :mip_suboptimal,
    XPRS_MIP_INFEAS           => :mip_infeasible,
    XPRS_MIP_OPTIMAL          => :mip_optimal,
    XPRS_MIP_UNBOUNDED        => :mip_unbounded
)

@xprs_int_attr get_mip_status_code XPRS_MIPSTATUS
@xprs_int_attr get_lp_status_code XPRS_LPSTATUS

get_mip_status(model::Model) = status_symbols_mip[get_mip_status_code(model)]::Symbol
get_lp_status(model::Model) = status_symbols_lp[get_lp_status_code(model)]::Symbol

@xprs_int_attr get_sol_count XPRS_MIPSOLS
@xprs_int_attr get_node_count XPRS_NODES

@xprs_int_attr get_simplex_iter_count XPRS_SIMPLEXITER
@xprs_int_attr get_barrier_iter_count XPRS_BARITER

# objective
@xprs_dbl_attr get_lp_objval XPRS_LPOBJVAL
@xprs_dbl_attr get_mip_objval XPRS_MIPBESTOBJVAL
@xprs_dbl_attr get_mip_objval_last XPRS_MIPOBJVAL

@xprs_dbl_attr get_bardualobj XPRS_BARDUALOBJ
@xprs_dbl_attr get_barprimalobj XPRS_BARPRIMALOBJ

@xprs_dbl_attr get_bestbound XPRS_BESTBOUND

@xprs_dbl_attr get_objrhs XPRS_OBJRHS

#iis
@xprs_int_attr num_iis XPRS_NUMIIS

function get_objval(model::Model)
    if num_intents(model)+num_sos(model) > 0
        return get_mip_objval(model)
    else
        return get_lp_objval(model)
    end
end

mutable struct OptimInfo
    status_lp::Symbol
    status_mip::Symbol

    #time ??

    sol_count::Int
    node_count::Int

    simplex_iter_count::Int
    barrier_iter_count::Int

end

function get_optiminfo(model::Model)
    OptimInfo(
        get_lp_status(model),
        get_mip_status(model),

        #get_runtime(model),

        get_sol_count(model),
        get_node_count(model),

        get_simplex_iter_count(model),
        get_barrier_iter_count(model)
    )
end

function show(io::IO, s::OptimInfo)
    println(io, "Xpress Optimization Info")
    println(io, "    lp status    = $(s.status_lp)")
    println(io, "    mip status   = $(s.status_mip)")
    #println(io, "    runtime  = $(s.runtime)")
    println(io, "    # solutions     = $(s.sol_count)")
    println(io, "    # nodes         = $(s.node_count)")

    println(io, "    # simplex iters = $(s.simplex_iter_count)")
    println(io, "    # barrier iters = $(s.barrier_iter_count)")
end

#################################################
#
#  solution query
#
#################################################

function get_complete_lp_solution(model::Model)
    cols = num_vars(model)
    rows = num_constrs(model)

    x = Vector{Float64}(undef, cols)
    slack = Vector{Float64}(undef, rows)
    dual = Vector{Float64}(undef, rows)
    red = Vector{Float64}(undef, cols)

    ret = @xprs_ccall(getlpsol, Cint,
        (Ptr{Nothing},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64}
            ), model.ptr_model, x, slack, dual, red)
    if ret != 0
        throw(XpressError(model))
    end
    return x, slack, dual, red
end

"""
    get_lp_solution(model::Model)

return a vector with primal variable solutions
"""
function get_lp_solution(model::Model)
    cols = num_vars(model)
    x = Vector{Float64}(undef, cols)
    get_lp_solution!(model, x)
    return x
end

"""
    get_lp_solution!(model::Model, x::Vector{Float64})

return a vector with primal variable solutions - inplace
"""
function get_lp_solution!(model::Model, x::Vector{Float64})
    _chklen(x, num_vars(model))
    ret = @xprs_ccall(getlpsol, Cint,
        (Ptr{Nothing},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64}
            ), model.ptr_model, x, C_NULL, C_NULL, C_NULL)
    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end

"""
    get_lp_slack(model::Model)

return a vector of slack values for rows (some srt of constraint primal solution)
"""
function get_lp_slack(model::Model)
    rows = num_constrs(model)
    slack = Vector{Float64}(undef, rows)
    get_lp_slack!(model, slack)
    return slack
end

"""
    get_lp_slack!(model::Model, slack::Vector{Float64})

return a vector of slack values for rows (some srt of constraint primal solution) - inplace
"""
function get_lp_slack!(model::Model, slack::Vector{Float64})
    _chklen(slack, num_constrs(model))
    ret = @xprs_ccall(getlpsol, Cint,
        (Ptr{Nothing},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64}
            ), model.ptr_model, C_NULL, slack, C_NULL, C_NULL)
    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end
function get_lp_slack_lin!(model::Model, place::Vector{Float64})
    if num_qconstrs(model) > 0
        _chklen(place, num_linconstrs(model))

        slack = zeros(num_constrs(model))
        get_lp_slack!(model, slack)

        lrows = get_lrows(model)
        _chklen(lrows, num_linconstrs(model))

        for i in eachindex(lrows)
            place[i] = slack[lrows[i]]
        end
    else
        get_lp_slack!(model, place)
    end

    return nothing
end

"""
    get_lp_dual(model::Model)

Return a vector of constraint dual solution
"""
function get_lp_dual(model::Model)
    rows = num_constrs(model)

    dual = Vector{Float64}(undef, rows)

    get_lp_dual!(model, dual)
    return dual
end

"""
    get_lp_dual!(model::Model, dual::Vector{Float64})

Return a vector of constraint dual solution - inplace
"""
function get_lp_dual!(model::Model, dual::Vector{Float64})
    rows = num_constrs(model)
    _chklen(dual, rows)
    ret = @xprs_ccall(getlpsol, Cint,
        (Ptr{Nothing},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64}
            ), model.ptr_model, C_NULL, C_NULL, dual, C_NULL)
    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end
function get_lp_dual_lin!(model::Model, place::Vector{Float64})
    if num_qconstrs(model) > 0
        _chklen(place, num_linconstrs(model))

        dual = zeros(num_constrs(model))
        get_lp_dual!(model, dual)

        lrows = get_lrows(model)
        _chklen(lrows, num_linconstrs(model))

        for i in eachindex(lrows)
            place[i] = dual[lrows[i]]
        end
    else
        get_lp_dual!(model, place)
    end

    return nothing
end

"""
    get_lp_reducedcost(model::Model)

Return a vector of variable reduced cost (dual solution)
"""
function get_lp_reducedcost(model::Model)
    cols = num_vars(model)

    red = Vector{Float64}(undef, cols)

    get_lp_reducedcost!(model, red)
    return red
end

"""
    get_lp_reducedcost!(model::Model, red::Vector{Float64})

Return a vector of variable reduced cost (dual solution) - inplace
"""
function get_lp_reducedcost!(model::Model, red::Vector{Float64})
    cols = num_vars(model)

    _chklen(red, cols)

    ret = @xprs_ccall(getlpsol, Cint,
        (Ptr{Nothing},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64}
            ), model.ptr_model, C_NULL, C_NULL, C_NULL, red)
    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end

## mip sol
"""
    get_mip_solution(model::Model)

Return a vector of variable primal solutions
"""
function get_mip_solution(model::Model)
    cols = num_vars(model)

    x = Vector{Float64}(undef, cols)

    get_mip_solution!(model, x)
    return x
end

"""
    get_mip_solution!(model::Model, x::Vector{Float64})

Return a vector of variable primal solutions - inplace
"""
function get_mip_solution!(model::Model, x::Vector{Float64})
    cols = num_vars(model)

    _chklen(x, cols)

    ret = @xprs_ccall(getmipsol, Cint,
        (Ptr{Nothing},
         Ptr{Float64},
         Ptr{Float64}
            ), model.ptr_model, x, C_NULL)
    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end

"""
    get_mip_slack(model::Model)

Return a vector of constraint primal solutions (slacks)
"""
function get_mip_slack(model::Model)
    rows = num_constrs(model)

    slack = Vector{Float64}(undef, rows)

    get_mip_slack!(model, slack)
    return slack
end

"""
    get_mip_slack!(model::Model, slack::Vector{Float64})

Return a vector of constraint primal solutions (slacks) - inplace
"""
function get_mip_slack!(model::Model, slack::Vector{Float64})
    rows = num_constrs(model)
    _chklen(slack, rows)
    ret = @xprs_ccall(getmipsol, Cint,
        (Ptr{Nothing},
         Ptr{Float64},
         Ptr{Float64}
            ), model.ptr_model, C_NULL, slack)
    if ret != 0
        throw(XpressError(model))
    end
    return slack
end

function get_mip_slack_lin!(model::Model, place::Vector{Float64})
    if num_qconstrs(model) > 0
        _chklen(place, num_linconstrs(model))

        slack = zeros(num_constrs(model))
        get_mip_slack!(model, slack)

        lrows = get_lrows(model)
        _chklen(lrows, num_linconstrs(model))

        for i in eachindex(lrows)
            place[i] = slack[lrows[i]]
        end
    else
        get_mip_slack!(model, place)
    end
end

"""
    get_solution!(model::Model, x::Vector{Float64})

Return a vector of variable primal solution - inplace
"""
function get_solution!(model::Model, x::Vector{Float64})
    if is_mip(model)
        return get_mip_solution!(model,x)
    else
        return get_lp_solution!(model,x)
    end
end

"""
    get_solution(model::Model)

Return a vector of variable primal solution
"""
function get_solution(model::Model)
    if is_mip(model)
        return get_mip_solution(model)
    else
        return get_lp_solution(model)
    end
end

"""
    get_slack!(model::Model, slack::Vector{Float64})

Return a vector of constraints slacks - inplace

    get_slack(model::Model)

Return a vector of constraints slacks
"""
function get_slack!(model::Model, slack::Vector{Float64})
    if is_mip(model)
        return get_mip_slack!(model, slack)
    else
        return get_lp_slack!(model, slack)
    end
end
function get_slack(model::Model)
    if is_mip(model)
        return get_mip_slack(model)
    else
        return get_lp_slack(model)
    end
end
function get_slack_lin!(model::Model, slack::Vector{Float64})
    if is_mip(model)
        return get_mip_slack_lin!(model, slack)
    else
        return get_lp_slack_lin!(model, slack)
    end
end

"""
    get_dual!(model::Model, dual::Vector{Float64})

Return a vector of constraint dual solutions
"""
function get_dual!(model::Model, dual::Vector{Float64})
    if is_mip(model)
        error("Not possible to get MIP duals")
        return Float64[]
    else
        return get_lp_dual!(model,dual)
    end
end
function get_dual_lin!(model::Model, dual::Vector{Float64})
    if is_mip(model)
        error("Not possible to get MIP duals")
        return Float64[]
    else
        return get_lp_dual_lin!(model,dual)
    end
end

"""
    get_dual(model::Model)

Return a vector of constraint dual solutions
"""
function get_dual(model::Model)
    if is_mip(model)
        error("Not possible to get MIP duals")
        return Float64[]
    else
        return get_lp_dual(model)
    end
end

"""
    get_reducedcost!(model::Model, red::Vector{Float64})

Return a vector of variable dual solution
"""
function get_reducedcost!(model::Model, red::Vector{Float64})
    if is_mip(model)
        error("Not possible to get MIP reduced costs")
        return Float64[]
    else
        return get_lp_reducedcost!(model, red)
    end
end

"""
    get_reducedcost(model::Model)

Return a vector of variable dual solution
"""
function get_reducedcost(model::Model)
    if is_mip(model)
        error("Not possible to get MIP reduced costs")
        return Float64[]
    else
        return get_lp_reducedcost(model)
    end
end

#################################################
#
#  basis query
#
#################################################
const basicmap = Dict(
    0 => :NonbasicAtLower,#_or_SuperbasicAtZero,
    1 => :Basic,
    2 => :NonbasicAtUpper,
    3 => :Superbasic
)
const basicmap_rev = Dict(
    :NonbasicAtLower => Cint(0),#_or_SuperbasicAtZero,
    :Basic => Cint(1),
    :NonbasicAtUpper => Cint(2),
    :Superbasic => Cint(3)
)

"""
    loadbasis(model::Model, x::Vector)

Load basis to a problem in form of primal solutions

    loadbasis(model::Model, rval::Vector{Symbol}, cval::Vector{Symbol})

Load basis to a problem in terms of basicness description
Variables and constraints ca be: `Basic`, `NonBasicAtLower`, `NonBasicAtUpper` and `Superbasic`
"""
function loadbasis(model::Model, x::Vector)#, status::Symbol = :unstarted, isnew::Vector{Bool} = [false])

    ncols = num_vars(model)
    nrows = num_constrs(model)

    length(x) != ncols && error("solution candidate size is different from the number of columns")

    cvals = Array{Cint}(undef,  ncols)
    rvals = Array{Cint}(undef,  nrows)

    # obtain situation of columns

    lb = get_lb(model)
    ub = get_ub(model)

    for i in 1:ncols
        if isapprox(x[i],lb[i])
            cvals[i] = basicmap_rev[:NonbasicAtLower]
        elseif isapprox(x[i],ub[i])
            cvals[i] = basicmap_rev[:NonbasicAtUpper]
        else
            cvals[i] = basicmap_rev[:Basic]
        end
    end

    # obtain situation of rows where: y = Ax

    A = get_constrmatrix(model) #A is sparse

    y = A*x

    senses = get_rowtype(model)
    rhs    = get_rhs(model)

    for j in 1:nrows
        if senses[j] == XPRS_EQ && isapprox(y[j], rhs[j])
            rvals[j] = basicmap_rev[:NonbasicAtLower]
        elseif senses[j] == XPRS_GEQ && isapprox(y[j], rhs[j])
            rvals[j] = basicmap_rev[:NonbasicAtLower]
        elseif senses[j] == XPRS_LEQ && isapprox(y[j], rhs[j])
            rvals[j] = basicmap_rev[:NonbasicAtUpper]
        else
            rvals[j] = basicmap_rev[:Basic]
        end
    end

    loadbasis(model, rvals, cvals)

    return nothing
end
function loadbasis(model::Model, rval::Vector{Symbol}, cval::Vector{Symbol})

    nrval = map(x->basicmap_rev[x], rval)
    ncval = map(x->basicmap_rev[x], cval)

    loadbasis(model, nrval, ncval)

    return nothing
end
function loadbasis(model::Model, rval::Vector{Cint}, cval::Vector{Cint})
# int XPRS_CC XPRSloadbasis(XPRSprob prob, const int rstatus[], const intcstatus[]);

    ret = @xprs_ccall(loadbasis, Cint,
        (Ptr{Nothing},
         Ptr{Cint},
         Ptr{Cint}
            ), model.ptr_model, rval, cval)
    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end

function get_basis(model::Model)
    cval = Array{Cint}(undef,  num_vars(model))
    cbasis = Array{Symbol}( num_vars(model))

    rval = Array{Cint}(undef,  num_constrs(model))
    rbasis = Array{Symbol}( num_constrs(model))

    ret = @xprs_ccall(getbasis, Cint,
        (Ptr{Nothing},
         Ptr{Cint},
         Ptr{Cint}
            ), model.ptr_model, rval, cval)
    if ret != 0
        throw(XpressError(model))
    end

    for it in 1:length(cval)
        cbasis[it] = basicmap[cval[it]]
    end

    for it in 1:length(rval)
        rbasis[it] = basicmap[rval[it]]
    end
    return cbasis, rbasis
end

#################################################
#
#  IIS query
#
#################################################

function get_iisdata(model::Model, num::Int)
# num is the number of THE IIS to be queried

    rows = Array{Cint}(undef, 1)
    cols = Array{Cint}(undef, 1)
    ret = @xprs_ccall(getiisdata, Cint,
        (Ptr{Nothing},
            Cint,# num
         Ptr{Cint},# #row
         Ptr{Cint},# #cols
         Ptr{Cint},# rows set
         Ptr{Cint},# bounds set
         Ptr{Cchar},# cstr type
         Ptr{Cchar},# bnd type
         Ptr{Float64}, # duals
         Ptr{Float64}, # reduced
         Ptr{Cchar}, # isol rows
         Ptr{Cchar} #isol cols
            ), model.ptr_model, Cint(num), rows,
        cols,C_NULL,C_NULL,
        C_NULL,C_NULL,C_NULL,
        C_NULL,C_NULL,C_NULL,)

    if ret != 0
        throw(XpressError(model))
    end

    rows_set = Array{Cint}(undef,  rows[1])
    cols_set = Array{Cint}(undef,  cols[1])

    ret = @xprs_ccall(getiisdata, Cint,
        (Ptr{Nothing},
            Cint,# num
         Ptr{Cint},# #row
         Ptr{Cint},# #cols
         Ptr{Cint},# rows set
         Ptr{Cint},# bounds set
         Ptr{Cchar},# cstr type
         Ptr{Cchar},# bnd type
         Ptr{Float64}, # duals
         Ptr{Float64}, # reduced
         Ptr{Cchar}, # isol rows
         Ptr{Cchar} #isol cols
            ), model.ptr_model, Cint(num), rows,
        cols,rows_set,cols_set,
        C_NULL,C_NULL,C_NULL,
        C_NULL,C_NULL,C_NULL,)

    if ret != 0
        throw(XpressError(model))
    end

    return rows_set + 1, cols_set + 1
end

#################################################
#
#  infeasibility & unboundedness
#
#################################################

function hasdualray(model::Model)::Bool

    hasray = Array{Cint}(undef,  1)

    ret = @xprs_ccall(getdualray, Cint,
        (Ptr{Nothing},
            Ptr{Float64},
            Ptr{Cint}
            ), model.ptr_model, C_NULL, hasray)
    if ret != 0
        throw(XpressError(model))
    end

    if hasray[1] == 0
        return false
    else
        return true
    end
end

function getdualray!(model::Model, ray::Vector{Float64})

    @assert length(ray) == num_constrs(model)

    hasray = Array{Cint}(undef, 1)

    ret = @xprs_ccall(getdualray, Cint,
        (Ptr{Nothing},
            Ptr{Float64},
            Ptr{Cint}
            ), model.ptr_model, ray, hasray)
    if ret != 0
        throw(XpressError(model))
    end
    return hasray[1] != 0
end

function getdualray(model::Model)

    dray = Array{Float64}(undef,  num_constrs(model))

    if !getdualray!(model, dray)
        Compat.@warn("Xpress solver was unable to provide an infeasibility ray")
        return dray
    end

    return dray
end

function hasprimalray(model::Model)::Bool

    hasray = Array{Cint}(undef,  1)

    ret = @xprs_ccall(getprimalray, Cint,
        (Ptr{Nothing},
            Ptr{Float64},
            Ptr{Cint}
            ), model.ptr_model, C_NULL, hasray)
    if ret != 0
        throw(XpressError(model))
    end

    if hasray[1] == 0
        return false
    else
        return true
    end
end

function getprimalray!(model::Model, ray::Vector{Float64})

    @assert length(ray) == num_vars(model)

    hasray = Array{Cint}(undef,  1)

    ret = @xprs_ccall(getprimalray, Cint,
        (Ptr{Nothing},
            Ptr{Float64},
            Ptr{Cint}
            ), model.ptr_model, ray, hasray)
    if ret != 0
        throw(XpressError(model))
    end
    return hasray[1] != 0
end

function getprimalray(model::Model)

    pray = Array{Float64}(undef,  num_vars(model))

    if !getprimalray!(model, pray)
        Compat.@warn("Xpress solver was unable to provide an unboundedness ray")
        return pray
    end

    return pray
end

function repairweightedinfeasibility(model::Model, scode::Vector{Cint}, lrp::Vector{Float64}, grp::Vector{Float64}, lbp::Vector{Float64}, ubp::Vector{Float64}, phase2::Cchar = Cchar('f'), delta::Float64=0.001, flags::String="")
# int XPRS_CC XPRSrepairweightedinfeas(XPRSprob prob, int *scode, const double lrp[], const double grp[], const double lbp[], const double ubp[], char phase2, double delta, const char *optflags)
    ret = @xprs_ccall(repairweightedinfeas, Cint,
        (Ptr{Nothing},
         Ptr{Cint},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64},
         Ptr{Float64},
         Cchar,
         Float64,
         Ptr{UInt8}),
        model.ptr_model, scode, lrp, grp, lbp, ubp, phase2, delta, flags)
    if ret != 0
        throw(XpressError(model))
    end

    return nothing
end

"""
    repairweightedinfeasibility(model::Model, lrp::Vector{Float64}, grp::Vector{Float64}, lbp::Vector{Float64}, ubp::Vector{Float64}; phase2::Cchar = Cchar('f'), delta::Float64=0.001, flags::String="")
    repairweightedinfeasibility(model::Model, phase2::Cchar = Cchar('f'), delta::Float64=0.001, flags::String="")

Repair infeasibility tool
"""
repairweightedinfeasibility(model::Model, lrp::Vector{Float64}, grp::Vector{Float64}, lbp::Vector{Float64}, ubp::Vector{Float64}; phase2::Cchar = Cchar('f'), delta::Float64=0.001, flags::String="") = repairweightedinfeasibility(model::Model, lrp::Vector{Float64}, grp::Vector{Float64}, lbp::Vector{Float64}, ubp::Vector{Float64}, phase2, delta, flags)

function repairweightedinfeasibility(model::Model, lrp::Vector{Float64}, grp::Vector{Float64}, lbp::Vector{Float64}, ubp::Vector{Float64}, phase2::Cchar, delta::Float64, flags::String)

    cols = num_vars(model)
    rows = num_constrs(model)

    _chklen(rows, lrp)
    _chklen(rows, grp)
    _chklen(cols, lbp)
    _chklen(cols, ubp)

    scode = zeros(Cint,1)

    repairweightedinfeasibility(model, scode, lrp, grp, lbp, ubp, phase2, delta, flags)

    return scode[1]
end

function repairweightedinfeasibility(model::Model, phase2::Cchar = Cchar('f'), delta::Float64=0.001, flags::String="")
    cols = num_vars(model)
    rows = num_constrs(model)
    lrp = ones(rows)
    grp = ones(rows)
    lbp = ones(cols)
    ubp = ones(cols)

    return repairweightedinfeasibility(model, lrp, grp, lbp, ubp, phase2, delta, flags)
end

@xprs_int_attr stopstatus XPRS_STOPSTATUS

@enum StopStatus StopNone StopTimeLimit StopControlC StopNodeLimit StopIterLimit StopMIPGap StopSolLimit StopUser StopUnknown

function get_stopstatus(m::Model)
    ss = stopstatus(m)

    if ss == XPRS_STOP_NONE
        return StopNone
    elseif ss == XPRS_STOP_TIMELIMIT
        return StopTimeLimit
    elseif ss == XPRS_MAXTIME
        return StopTimeLimit
    elseif ss == XPRS_STOP_CTRLC
        return StopControlC
    elseif ss == XPRS_STOP_NODELIMIT
        return StopNodeLimit
    elseif ss == XPRS_STOP_ITERLIMIT
        return StopIterLimit
    elseif ss == XPRS_STOP_MIPGAP
        return StopMIPGap
    elseif ss == XPRS_STOP_SOLLIMIT
        return StopSolLimit
    elseif ss == XPRS_STOP_USER
        return StopUser
    end

    return StopUnknown
end


text:
using CSV, DataFrames, DataFramesMeta
respath = "benchmark/bench_string_results/"
paths = readdir(respath)
@time df = reduce(vcat, CSV.read.(respath.*paths))

df1 = by(df, [:test,:n,:strlen, :strlentype]) do df1
    mean(df1[:elapsed])
end
rename!(df1,:x1=>:elapsed)

df_var = df1[df1[:strlentype] .== "variable",:]
df_fix = df1[df1[:strlentype] .== "fixed",:]
unstack(df_var, [:n, :strlen],:test,:elapsed)
unstack(df_fix, [:n, :strlen],:test,:elapsed)

df_id = df1[df1[:strlentype] .== "id",:]
unstack(df_id, [:n, :strlen],:test,:elapsed)

# find the winnder of fixed length string
function findwinnerstrsort(df1)
    df_winner = by(df1,[:n,:strlen]) do subdf
        me = minimum(subdf[:elapsed])
        subdf[subdf[:elapsed] .== me,[:test,:elapsed]]
    end

    df_no_r = df1[df1[:test] .!= "r_radixsort_elapsed",:]
    df_winner_no_r = by(df_no_r,[:n,:strlen]) do subdf
        me = minimum(subdf[:elapsed])
        subdf[subdf[:elapsed] .== me,[:test,:elapsed]]
    end

    rename!(df_winner_no_r,:test=>:test2)
    rename!(df_winner_no_r,:elapsed=>:elapsed2)

    dfres = join(df_winner, df_winner_no_r, on=[:n,:strlen])
    dfres[:speed_ratio] = dfres[:elapsed2]./dfres[:elapsed]
    dfres
end

findwinnerstrsort(df_var)
findwinnerstrsort(df_fix)
findwinnerstrsort(df_id)

ext: jl


text:
export iauFapa03
"""
Fundamental argument, IERS Conventions (2003):
general accumulated precession in longitude.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   t     double    TDB, Julian centuries since J2000.0 (Note 1)

Returned (function value):
         double    general precession in longitude, radians (Note 2)

Notes:

   1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

   2. The expression used is as adopted in IERS Conventions (2003).  It
   is taken from Kinoshita & Souchay (1990) and comes originally
   from Lieske et al. (1977).

References:

   Kinoshita, H. and Souchay J. 1990, Celest.Mech. and Dyn.Astron.
   48, 187

   Lieske, J.H., Lederle, T., Fricke, W. & Morando, B. 1977,
   Astron.Astrophys. 58, 1-16

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauFapa03(double t)
function iauFapa03(t::Real)
   return ccall((:iauFapa03, libsofa_c), Cdouble, (Cdouble,), convert(Float64, t))
end


text:
# export HJICache, GriddedInterpolation, Gridded, Linear, SVector    # TODO: I shouldn't have to do this

@maintain_type struct SimpleCarState{T} <: FieldVector{4,T}
    E::T    # world frame "x" position of CM
    N::T    # world frame "y" position of CM
    ψ::T    # world frame heading of vehicle
    V::T    # speed of vehicle (velocity assumed to be in the heading direction)
end

@maintain_type struct HJIRelativeState{T} <: FieldVector{7,T}
    ΔE::T
    ΔN::T
    Δψ::T
    Ux::T
    Uy::T
    V::T
    r::T
end

function HJIRelativeState(us::BicycleState, them::SimpleCarState)
    cψ, sψ = sincos(us.ψ)
    ΔE, ΔN = @SMatrix([cψ sψ; -sψ cψ])*SVector(them.E - us.E, them.N - us.N)
    HJIRelativeState(ΔE, ΔN, adiff(them.ψ, us.ψ), us.Ux, us.Uy, them.V, us.r)
end

struct HJICache
    grid_knots::NTuple{7,Vector{Float32}}
    V ::GriddedInterpolation{Float32,7,Float32,Gridded{Linear},NTuple{7,Vector{Float32}}}
    ∇V::GriddedInterpolation{SVector{7,Float32},7,SVector{7,Float32},Gridded{Linear},NTuple{7,Vector{Float32}}}
end

function placeholder_HJICache()
    grid_knots = tuple((Float32[-1000., 1000.] for i in 1:7)...)
    V  = interpolate(Float32, Float32, grid_knots, zeros(Float32,2,2,2,2,2,2,2), Gridded(Linear()))
    ∇V = interpolate(Float32, SVector{7,Float32}, grid_knots, zeros(SVector{7,Float32},2,2,2,2,2,2,2), Gridded(Linear()))
    HJICache(grid_knots, V, ∇V)
end

function HJICache(fname::String)
    # if endswith(fname, ".mat")
    #     HJIdata = matread(fname)
    #     grid_knots = tuple((x -> convert(Vector{Float32}, vec(x))).(HJIdata["avoid_set"]["g"]["vs"])...)
    #     V  = interpolate(Float32, Float32, grid_knots,
    #                      HJIdata["avoid_set"]["data"][1,1], Gridded(Linear()))
    #     ∇V = interpolate(Float32, SVector{7,Float32}, grid_knots,
    #                      SVector.(HJIdata["avoid_set"]["deriv"]...), Gridded(Linear()))
    #     HJICache(grid_knots, V, ∇V)
    # else
    if endswith(fname, ".jld2")
        @load fname grid_knots V_raw ∇V_raw
        V  = interpolate(Float32, Float32, grid_knots, V_raw, Gridded(Linear()))
        ∇V = interpolate(Float32, SVector{7,Float32}, grid_knots, Array(reinterpret(SVector{7,Float32}, ∇V_raw)), Gridded(Linear()))
        HJICache(grid_knots, V, ∇V)
    else
        error("Unknown file type for loading HJICache")
    end
end

function save(fname::String, cache::HJICache)
    grid_knots = cache.grid_knots
    V_raw  = cache.V.coefs
    ∇V_raw = Array(reinterpret(Float32, cache.∇V.coefs))
    @save fname grid_knots V_raw ∇V_raw
end

function Base.getindex(cache::HJICache, x::HJIRelativeState{T}) where {T}
    if all(cache.grid_knots[i][1] <= x[i] <= cache.grid_knots[i][end] for i in 1:length(cache.grid_knots))
        (V=cache.V(x[1], x[2], x[3], x[4], x[5], x[6], x[7]), ∇V=cache.∇V(x[1], x[2], x[3], x[4], x[5], x[6], x[7]))    # avoid splatting penalty
    else
        (V=T(Inf), ∇V=zeros(SVector{7,T}))
    end
end

function relative_dynamics(X::VehicleModel, (ΔE, ΔN, Δψ, Ux, Uy, V, r)::StaticVector{7},    # relative state
                                            uR::StaticVector{2},                            # robot control
                                            (ω, a)::StaticVector{2})                        # human control
    bicycle_dynamics = X(BicycleState(ΔE, ΔN, Δψ, Ux, Uy, r), uR)
    sΔψ, cΔψ = sincos(Δψ)
    SVector(
        V*cΔψ - Ux + ΔN*r,
        V*sΔψ - Uy - ΔE*r,
        ω - r,
        bicycle_dynamics[4],
        bicycle_dynamics[5],
        a,
        bicycle_dynamics[6]
    )
end

function optimal_disturbance(X::VehicleModel, relative_state::HJIRelativeState, ∇V::StaticVector{7}, dMode=:min)
    BM, CL = X.bicycle_model, X.control_limits
    μ, m, G = BM.μ, BM.m, BM.G
    κ_max   = CL.κ_max
    Ax_max  = CL.Fx_max/m
    Pmx_max = CL.Px_max/m
    maxA = 0.9*μ*G

    sgn = (dMode == :max ? 1 : -1)

    V = relative_state.V
    lam_w  = ∇V[3]    # psi_rel dot
    lam_Ax = ∇V[6]    # v_dot
    lam_Ay = lam_w/V
    lam_norm = hypot(lam_Ax, lam_Ay)

    if lam_norm < 1e-3
        return SVector{2,Float64}(0, 0)
    else
        desAx = sgn*lam_Ax*maxA/lam_norm
        desAy = sgn*lam_Ay*maxA/lam_norm
        maxAx = min(Ax_max, Pmx_max/V)
        maxAy = κ_max*V*V
        if desAx > maxAx
            if abs(desAy) < maxAy
                maxAy = min(maxAy, sqrt(maxA*maxA - maxAx*maxAx))
            end
            return SVector{2,Float64}(copysign(maxAy, desAy)/V, maxAx)
        else
            if abs(desAy) > maxAy
                if desAx > 0
                    maxAx = min(sqrt(maxA*maxA - maxAy*maxAy), maxAx)
                    return SVector{2,Float64}(copysign(maxAy, desAy)/V, maxAx)
                else
                    return SVector{2,Float64}(copysign(maxAy, desAy)/V, -sqrt(maxA*maxA - maxAy*maxAy))
                end
            else
                return SVector{2,Float64}(desAy/V, maxAx)
            end
        end
    end
end

function optimal_control(X::VehicleModel, relative_state::HJIRelativeState, ∇V::StaticVector{7}, uMode=:max; N=50)
    BM, LP, CL = X.bicycle_model, X.longitudinal_params, X.control_limits
    L, a, b, h, m, μ, Cαf, Cαr, G, Izz = BM.L, BM.a, BM.b, BM.h, BM.m, BM.μ, BM.Cαf, BM.Cαr, BM.G, BM.Izz
    δ_max, Fx_max, Fx_min = CL.δ_max, CL.Fx_max, CL.Fx_min
    fake_qR = SVector(0., 0., 0., relative_state.Ux, relative_state.Uy, relative_state.r)    # sufficient for lateral_tire_forces

    sgn = (uMode == :max ? 1 : -1)
    A = ∇V[4]/m
    B = ∇V[5]/m + a*∇V[7]/Izz
    C = ∇V[5]/m - b*∇V[7]/Izz
    δ_opt  = ifelse(B >= 0, sgn*δ_max, -sgn*δ_max)
    V_opt  = -sgn*Inf
    Fx_opt = 0.0
    for n in 0:(N-1)
        frac = n/(N-1)
        Fx = frac*Fx_max + (1 - frac)*Fx_min
        uR = BicycleControl(δ_opt, longitudinal_tire_forces(LP, Fx)...)
        Fyf, Fyr = lateral_tire_forces(BM, fake_qR, uR)
        V = A*Fx + B*Fyf + C*Fyr
        if (sgn > 0 && V > V_opt) || (sgn < 0 && V < V_opt)
            Fx_opt = Fx
            V_opt  = V
        end
    end
    BicycleControl2(δ_opt, Fx_opt)
end

function compute_reachability_constraint(X::VehicleModel, cache::HJICache, relative_state::HJIRelativeState, ϵ,
                                         uR_lin=optimal_control(X, relative_state, cache[relative_state].∇V))    # definitely not the correct choice...
    V, ∇V = cache[relative_state]
    if V > ϵ
        (M=SVector{2,Float64}(0, 0), b=1.0)
    else
        uH_opt = optimal_disturbance(X, relative_state, ∇V)
        ∇H_uR = ForwardDiff.gradient(uR -> dot(∇V, relative_dynamics(X, SVector(relative_state), uR, SVector(uH_opt))), SVector(uR_lin))
        (M=∇H_uR, b=dot(∇V, relative_dynamics(X, relative_state, uR_lin, uH_opt)) - dot(∇H_uR, uR_lin)) # so that H = dot(∇V, uR) ≈ ∇H_uR*uR + c
    end
end

max_stars_repo_name: Michelle-NYX/CollisionAvoidance

max_issues_repo_name: Michelle-NYX/CollisionAvoidance

max_forks_repo_name: Michelle-NYX/CollisionAvoidance


text:
struct SDP <: AbstractLoss end

function complete(D::MaskedEDM{T}, loss::SDP, dims::Int) where T <: Real
    return SDP(D)
    # return EuclideanDistanceMatrix(D, distances=true)
end

function SDP(D::MaskedEDM{T}) where T <: Real
    ##
    # D = sdr_complete_D(D, W, dims)
    #
    # Completes a partially observed D (with observed entries assumed noiseless)
    # using the semidefinite relaxation (SDR).
    #
    # INPUT:  D ... input partial EDM
    #         W   ... observation mask (1 for observed, 0 for non-observed entries)
    #
    # OUTPUT: D   ... completed EDM
    #
    # Author: Ivan Dokmanic, 2014

    Dm = D.D .* D.mask
    n = size(Dm, 1)

    # The old SVD code to do the following is also good.
    x = -1/(n + sqrt(n))
    y = -1/sqrt(n)
    V = [y * ones(1, n-1); x * ones(n-1, n-1) + Matrix(I, n-1, n-1)]
    e = ones(n, 1)

    G = Semidefinite(n-1)
    B = V * G * V'
    E = diag(B) * e' + e * diag(B)' - 2 * B

    problem = maximize(tr(G))
    problem.constraints += E .* D.mask == D.D .* D.mask
    problem.constraints += isposdef(G)

    solve!(problem, SCSSolver(verbose=false))

    # Can I just use what's above inside the cvx block? Don't do that, do rank
    # thresholding here!
    B = V * G.value * V'
    D = B[diagind(B)] * e' + e * B[diagind(B)]' .- 2*B

    # return EuclideanDistanceMatrix(D, distances=true)
    return D
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/nghttp2_jll.jl
module nghttp2_jll

using Libdl

const PATH_list = String[]
const LIBPATH_list = String[]

export libnghttp2

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libnghttp2_handle = C_NULL
libnghttp2_path = ""

if Sys.iswindows()
    const libnghttp2 = "libnghttp2-14.dll"
elseif Sys.isapple()
    const libnghttp2 = "@rpath/libnghttp2.14.dylib"
else
    const libnghttp2 = "libnghttp2.so.14"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libnghttp2_handle = dlopen(libnghttp2)
    global libnghttp2_path = dlpath(libnghttp2_handle)
end

is_available() = true

end  # module nghttp2_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
using ParametricInversion
using Pkg
Pkg.develop(path = joinpath(pwd(), "ParametricInversionTestUtils"))

using ParametricInversionTestUtils
using Test
using Random
Random.seed!(0)

@testset "ParametricInversion.jl" begin
  include("invert.jl")
end


text:
#!/usr/bin/env julia

export pe8

module PE8

const data=parse.(Int,collect(
    "73167176531330624919225119674426574742355349194934"*
    "96983520312774506326239578318016984801869478851843"*
    "85861560789112949495459501737958331952853208805511"*
    "12540698747158523863050715693290963295227443043557"*
    "66896648950445244523161731856403098711121722383113"*
    "62229893423380308135336276614282806444486645238749"*
    "30358907296290491560440772390713810515859307960866"*
    "70172427121883998797908792274921901699720888093776"*
    "65727333001053367881220235421809751254540594752243"*
    "52584907711670556013604839586446706324415722155397"*
    "53697817977846174064955149290862569321978468622482"*
    "83972241375657056057490261407972968652414535100474"*
    "82166370484403199890008895243450658541227588666881"*
    "16427171479924442928230863465674813919123162824586"*
    "17866458359124566529476545682848912883142607690042"*
    "24219022671055626321111109370544217506941658960408"*
    "07198403850962455444362981230987879927244284909188"*
    "84580156166097919133875499200524063689912560717606"*
    "05886116467109405077541002256983155200055935729725"*
    "71636269561882670428252483600823257530420752963450"
));

end

pe8(data::Vector{<:Integer}=PE8.data,n::Integer=13)=maximum(prod(data[i:i+n-1]) for i in 1:(length(data)-(n-1)))

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe8()|>print
end

max_stars_repo_path: 8. Largest product in a series/8.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 8. Largest product in a series/8.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 8. Largest product in a series/8.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
module FreudenthalTriangulation

using LinearAlgebra
using SparseArrays

export freudenthal_vertices, freudenthal_simplex, barycentric_coordinates, freudenthal_simplex_and_coords,
        freudenthal_simplex_and_coords!

# Construct the Freudenthal triangulation of the belief simplex

"""
    freudenthal_vertices(n::Int64, m::Int64)
Construct the list of Freudenthal vertices in an `n` dimensional space with grid resolution `m`.
The vertices are represented by a list of `n` dimensional vectors.
"""
function freudenthal_vertices(n::Int64, m::Int64)
    V = Vector{Int}[]
    v = Vector{Int}(undef, n)
    v[1] = m
    freudenthal_vertices!(V, v, 2)
    return V
end

function freudenthal_vertices!(V::Vector{Vector{Int64}}, v::Vector{Int64}, i::Int64)
    n = length(v)
    if i > n
        push!(V, copy(v))
        return
    end
    for k in 0 : v[i-1]
        v[i] = k
        freudenthal_vertices!(V, v, i+1)
    end
end

"""
    freudenthal_simplex(x::Vector{Float64})
Returns the list of vertices of the simplex of point `x` in the Freudenthal grid.
"""
function freudenthal_simplex(x::Vector{Float64})
    n = length(x)
    V = Vector{Vector{Int}}(undef, n+1)
    V[1] = floor.(Int, x)
    d = x - V[1]
    p = sortperm(d, rev=true)
    for i in 2 : n+1
        V[i] = copy(V[i-1])
        V[i][p[i-1]] += 1
    end
    return V
end

"""
    barycentric_coordinates(x::Vector{Int64}, V::Vector{Vector{Int64}})
Given a point `x` and its simplex `V` in the Freudenthal grid, returns the barycentric coordinates
of `x` in the grid. `V` must be in the same order as provided by the output of `freudenthal_simplex`
"""
function barycentric_coordinates(x, V)
    d = x - V[1]
    p = sortperm(d, rev=true)
    n = length(x)
    λ = Vector{Float64}(undef, n+1)
    λ[n+1] = d[p[n]]
    for i in n:-1:2
        λ[i] = d[p[i-1]] - d[p[i]]
    end
    λ[1] = 1.0 - sum(λ[2:end])
    return λ
end

"""
    freudenthal_coords(x::AbstractArray{Float64})
Given a point `x`, returns the simplex of the point `x` and the barycentric coordinates of `x` in the grid.
"""
function freudenthal_simplex_and_coords(x::Vector{Float64})
    V = freudenthal_simplex(x)
    return V, barycentric_coordinates(x, V)
end

"""
    freudenthal_simplex_and_coords!(x::AbstractArray{Float64}, V::Vector{Vector{Int64}}, λ::Vector{Float64})
Fills `V` and `λ` with the simplex points in the Freudenthal space and associated coordinates respectively.
"""
function freudenthal_simplex_and_coords!(x::AbstractArray{Float64}, V::Vector{Vector{Int64}}, λ::Vector{Float64})
    n = length(x)
    V[1] = floor.(Int, x)
    d = x - V[1]
    p = sortperm(d, rev=true)
    for i in 2 : n+1
        copyto!(V[i], V[i-1])
        V[i][p[i-1]] += 1
    end
    λ[n+1] = d[p[n]]
    for i in n:-1:2
        λ[i] = d[p[i-1]] - d[p[i]]
    end
    λ[1] = 1.0 - sum(λ[2:end])
    return V, λ
end

end # module


text:
mutable struct LogNormalFwdRateEuler{M <: AbstractMarketModel, G <: BrownianGenerator} <: AbstractMarketModelEvolver
  marketModel::M
  numeraires::Vector{Int}
  initialStep::Int
  generator::G
  fixedDrifts::Vector{Vector{Float64}}
  numberOfRates::Int
  numberOfFactors::Int
  curveState::LMMCurveState
  currentStep::Int
  forwards::Vector{Float64}
  displacement::Vector{Float64}
  logForwards::Vector{Float64}
  initialLogForwards::Vector{Float64}
  drifts1::Vector{Float64}
  initialDrifts::Vector{Float64}
  brownians::Vector{Float64}
  correlatedBrownians::Vector{Float64}
  alive::Vector{Int}
  calculators::Vector{LMMDriftCalculator}
end

function LogNormalFwdRateEuler(marketModel::AbstractMarketModel, factory::BrownianGeneratorFactory, numeraires::Vector{Int}, initialStep::Int = 1)
  numberOfRates = marketModel.numberOfRates
  numberOfFactors = marketModel.numberOfFactors
  curveState = LMMCurveState(marketModel.evolution.rateTimes)
  forwards = copy(marketModel.initialRates)
  displacements = copy(marketModel.displacements)
  logForwards = Vector{Float64}(undef, numberOfRates)
  initialLogForwards = Vector{Float64}(undef, numberOfRates)
  drifts1 = Vector{Float64}(undef, numberOfRates)
  initialDrifts = Vector{Float64}(undef, numberOfRates)
  brownians = Vector{Float64}(undef, numberOfFactors)
  correlatedBrownians = Vector{Float64}(undef, numberOfRates)
  alive = marketModel.evolution.firstAliveRate

  check_compatibility(marketModel.evolution, numeraires)

  steps = number_of_steps(marketModel.evolution)

  generator = create(factory, numberOfFactors, steps - (initialStep-1))

  currentStep = initialStep

  calculators = Vector{LMMDriftCalculator}(undef, steps)
  fixedDrifts = Vector{Vector{Float64}}(undef, steps)

  @inbounds @simd for j = 1:steps
    A = marketModel.pseudoRoots[j]
    calculators[j] = LMMDriftCalculator(A, displacements, marketModel.evolution.rateTaus, numeraires[j], alive[j])
    fixed = Vector{Float64}(undef, numberOfRates)

    for k in eachindex(fixed)
      # variance_ = dot(A[:, k], A[:, k])
      variance_ = dot(A[k, :], A[k, :])
      fixed[k] = -0.5 * variance_
    end

    fixedDrifts[j] = fixed
  end

  lognorm = LogNormalFwdRateEuler(marketModel, numeraires, initialStep, generator, fixedDrifts, numberOfRates, numberOfFactors, curveState, currentStep, forwards,
                              displacements, logForwards, initialLogForwards, drifts1, initialDrifts, brownians, correlatedBrownians, alive, calculators)

  set_forwards!(lognorm, marketModel.initialRates)
  return lognorm
end

function set_forwards!(lognorm::LogNormalFwdRateEuler, forwards::Vector{Float64})
  length(forwards) == lognorm.numberOfRates || error("mismatch between forwards and rate times")
  @inbounds @simd for i in eachindex(lognorm.initialLogForwards)
    lognorm.initialLogForwards[i] = log(forwards[i] + lognorm.displacement[i])
  end

  compute!(lognorm.calculators[lognorm.initialStep], forwards, lognorm.initialDrifts)

  return lognorm
end

function start_new_path!(lognorm::LogNormalFwdRateEuler)
  lognorm.currentStep = lognorm.initialStep
  lognorm.logForwards = copy(lognorm.initialLogForwards)

  return next_path!(lognorm.generator)
end

function advance_step!(lognorm::LogNormalFwdRateEuler)
  # we are going from T1 to T2

  # a) compute drifts D1 at T1
  if lognorm.currentStep > lognorm.initialStep
    compute!(lognorm.calculators[lognorm.currentStep], lognorm.forwards, lognorm.drifts1)
  else
    lognorm.drifts1 = copy(lognorm.initialDrifts)
  end

  # b) evolve forwards up to T2 using D1
  weight = next_step!(lognorm.generator, lognorm.brownians)
  A = copy(lognorm.marketModel.pseudoRoots[lognorm.currentStep])
  fixedDrift = copy(lognorm.fixedDrifts[lognorm.currentStep])

  alive = lognorm.alive[lognorm.currentStep]
  @inbounds @simd for i = alive:lognorm.numberOfRates
    lognorm.logForwards[i] += lognorm.drifts1[i] + fixedDrift[i]
    lognorm.logForwards[i] += dot(A[i, :], lognorm.brownians)
    lognorm.forwards[i] = exp(lognorm.logForwards[i]) - lognorm.displacement[i]
  end

  # same as PC evolver with two steps dropped

  # c) update curve state
  set_on_forward_rates!(lognorm.curveState, lognorm.forwards)
  lognorm.currentStep += 1

  return weight
end

current_state(lognorm::LogNormalFwdRateEuler) = lognorm.curveState
brownians_this_step(lognorm::LogNormalFwdRateEuler) = lognorm.brownians

# Clone #
clone(lognorm::LogNormalFwdRateEuler) = LogNormalFwdRateEuler(clone(lognorm.marketModel), copy(lognorm.numeraires), lognorm.initialStep, clone(lognorm.generator),
                                      deepcopy(lognorm.fixedDrifts), lognorm.numberOfRates, lognorm.numberOfFactors, clone(lognorm.curveState),
                                      lognorm.currentStep, copy(lognorm.forwards), copy(lognorm.displacement), copy(lognorm.logForwards), copy(lognorm.initialLogForwards),
                                      copy(lognorm.drifts1), copy(lognorm.initialDrifts), copy(lognorm.brownians), copy(lognorm.correlatedBrownians),
                                      copy(lognorm.alive), deepcopy(lognorm.calculators))


text:
@testset "Lazy derivatives" begin
    for (t₀,ϕ) in [(1.0,0.0), (4.0,π/3)]
        B = FEDVR(1.0:7, 4, t₀=t₀, ϕ=ϕ)
        B̃ = B[:,2:end-1]
        D = Derivative(axes(B,1))

        BD = B'D
        BDD = B'D'D

        @test BD*B == B'*D*B
        @test BDD*B == B'*D'*D*B

        @test B'D*B == B'*D*B
        @test B'D'D*B == B'*D'*D*B
    end
end

@testset "Materialize derivatives" begin
    @testset "$style materialization" for (style,first_derivative,second_derivative) in [
        ("Infix", (B,D) -> B'*D*B, (B,D) -> B'*D'*D*B),
        ("applied", (B,D) -> materialize(applied(*, B', D, B)),
         (B,D) -> materialize(applied(*, B', D', D, B))),
        ("apply", (B,D) -> apply(*, B', D, B),
         (B,D) -> apply(*, B', D', D, B))
    ]
        B = FEDVR(1.0:7, 4)
        D = Derivative(axes(B,1))

        ∇ = first_derivative(B,D)
        ∇² = second_derivative(B,D)

        A′ = Matrix(undef, B)
        A′′ = Matrix(undef, B)

        CompactBases.derop!(A′, B, 1)
        CompactBases.derop!(A′′, B, 2)

        @test ∇ isa BlockSkylineMatrix
        @test ∇² isa BlockSkylineMatrix

        @test ∇ == A′
        @test ∇² == A′′
    end
end

@testset "Derivatives in restricted bases" begin
    @testset "$style materialization" for (style,first_derivative,second_derivative) in [
        ("Infix", (B,D) -> B'*D*B, (B,D) -> B'*D'*D*B),
        ("applied", (B,D) -> materialize(applied(*, B', D, B)),
         (B,D) -> materialize(applied(*, B', D', D, B))),
        ("apply", (B,D) -> apply(*, B', D, B),
         (B,D) -> apply(*, B', D', D, B))
    ]
        @testset "order = $(order), sel = $(sel), N = $N" for (order,sel,N) in [
            (4,2:12,5),(4,1:13,5),([5,4,2,2],5:9,5),([5,4,3,2],5:10,5),
            (2,2:3,5),(2,2:4,5),(2,2:5,5),(2,1:5,5),
            (5,2:4,2)
        ]
            t = range(0,stop=20,length=N)

            R = FEDVR(t, order)
            R̃ = R[:,sel]

            D = Derivative(axes(R̃,1))

            expT = (any(order .> 2) ? BlockSkylineMatrix : Tridiagonal)

            ∇ = first_derivative(R, D)
            ∇̃ = first_derivative(R̃, D)

            @test ∇ isa expT
            @test ∇̃ isa expT

            @test Matrix(∇)[sel,sel] == Matrix(∇̃)

            ∇² = second_derivative(R, D)
            ∇̃² = second_derivative(R̃, D)

            @test ∇² isa expT
            @test ∇̃² isa expT

            @test Matrix(∇²)[sel,sel] == Matrix(∇̃²)
        end
    end

    @testset "Gradually decreasing size" begin
        L = 20.0
        n = 10
        t = range(-L/2, stop=L/2, length=n)
        B = FEDVR(t, [4,3,4,2,5,4,2,4,8])
        x = axes(B,1)

        D = Derivative(x)
        Dm = apply(*, B', D, B)
        @test Dm isa BlockSkylineMatrix

        @testset "From the left, a = $a" for a = 1:size(B,2)
            B̃ = B[:,a:end]
            D̃m = apply(*, B̃', D, B̃)
            @test D̃m isa BlockSkylineMatrix
            @test Dm[a:end,a:end] == D̃m
        end

        @testset "From the right, b = $b" for b = size(B,2):-1:1
            B̃ = B[:,1:b]
            D̃m = apply(*, B̃', D, B̃)
            @test D̃m isa BlockSkylineMatrix
            @test Dm[1:b,1:b] == D̃m
        end
    end
end

function compute_fedvr_derivative_errors(a, b, Ns, order::Integer, s::Integer, e::Integer,
                                         f::Function, g::Function, h::Function,
                                         verbosity=0)
    compute_derivative_errors(Ns, f, g, h, verbosity) do N
        t = range(a, stop=b, length=N)
        R = FEDVR(t, order)
        if (s,e) != (0,0)
            n = size(R,2)
            R = R[:,(1+s):(n-e)]
        end
        R
    end
end

@testset "Derivative accuracy" begin
    (f,g,h,a,b),s,e = derivative_test_functions(1.0), 1, 1

    orders = [(2,1.5,3),
              (3,1.2,2.5),
              (4,1.3,2.5),
              (5,1.3,2.5),
              (6,1.3,2.5),
              (7,1,2.5),
              (8,1,2.5),
              (9,1,2.5),
              (10,1,2.5)]
    slopes = zeros(length(orders),3)

    @info "FE-DVR derivative accuracy"
    for (i,(order,oa,ob)) in enumerate(orders)
        @show order
        Ns = ceil.(Int, 10 .^ range(oa,stop=ob,length=10))
        @show extrema(Ns)
        hs,ϵg,ϵh,ϵh′,pg,ph,ph′ = compute_fedvr_derivative_errors(a, b, Ns, order, s, e, f, g, h, 1)
        slopes[i,:] = [pg ph ph′]
    end

    oo = first.(orders)

    println("Derivative convergence rates:")
    pretty_table([oo slopes], ["Order", "pg", "ph", "ph′"])
    println()

    plt = lineplot(oo, slopes[:,1], name="pg", xlabel="FE-DVR order", ylabel="Error slope")
    lineplot!(plt, oo, slopes[:,2], name="ph")
    lineplot!(plt, oo, slopes[:,3], name="ph′")
    display(plt)
    println()
    println()

    # The convergence rate should be order - 1 (polynomial order =
    # polynomial degree + 1), but since the error slope fitting is a
    # bit error prone, we require that it is greater than order - 2.5.
    @test all(slopes[:,1] .> oo .- 2.5)
    @test all(slopes[:,1] .> oo .- 3.5)
    # Since the approximation to h is calculated by computing ∇∇f, we
    # lose one order extra, compared to ∇²f.
    @test all(slopes[:,2] .> oo .- 3.5)
    @test all(slopes[:,2] .> oo .- 4.5)
    @test all(slopes[:,3] .> oo .- 2.5)
    @test all(slopes[:,3] .> oo .- 4.0)
end


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

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

"""

Checks whether the specified account name is available or taken.
Param: subscriptionId::String (required)
Param: location::String (required)
Param: parameters::CheckNameAvailabilityParameters (required)
Param: api_version::String (required)
Return: NameAvailabilityInformation
"""
function _swaggerinternal_accountsCheckNameAvailability(_api::AccountsApi, subscriptionId::String, location::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", NameAvailabilityInformation, "/subscriptions/{subscriptionId}/providers/Microsoft.DataLakeStore/locations/{location}/checkNameAvailability", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function accountsCheckNameAvailability(_api::AccountsApi, subscriptionId::String, location::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsCheckNameAvailability(_api, subscriptionId, location, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function accountsCheckNameAvailability(_api::AccountsApi, response_stream::Channel, subscriptionId::String, location::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsCheckNameAvailability(_api, subscriptionId, location, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Creates the specified Data Lake Store account.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: parameters::CreateDataLakeStoreAccountParameters (required)
Param: api_version::String (required)
Return: DataLakeStoreAccount
"""
function _swaggerinternal_accountsCreate(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", DataLakeStoreAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function accountsCreate(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsCreate(_api, subscriptionId, resourceGroupName, accountName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function accountsCreate(_api::AccountsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsCreate(_api, subscriptionId, resourceGroupName, accountName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified Data Lake Store account.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_accountsDelete(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function accountsDelete(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsDelete(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function accountsDelete(_api::AccountsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsDelete(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Attempts to enable a user managed Key Vault for encryption of the specified Data Lake Store account.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_accountsEnableKeyVault(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/enableKeyVault", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function accountsEnableKeyVault(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsEnableKeyVault(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function accountsEnableKeyVault(_api::AccountsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsEnableKeyVault(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified Data Lake Store account.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Return: DataLakeStoreAccount
"""
function _swaggerinternal_accountsGet(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DataLakeStoreAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function accountsGet(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsGet(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function accountsGet(_api::AccountsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsGet(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists the Data Lake Store accounts within the subscription. The response includes a link to the next page of results, if any.
Param: subscriptionId::String (required)
Param: api_version::String (required)
Param: filter::String
Param: top::Int32
Param: skip::Int32
Param: select::String
Param: orderby::String
Param: count::Bool
Return: DataLakeStoreAccountListResult
"""
function _swaggerinternal_accountsList(_api::AccountsApi, subscriptionId::String, api_version::String; filter=nothing, top=nothing, skip=nothing, select=nothing, orderby=nothing, count=nothing, _mediaType=nothing)
    Swagger.validate_param("top", "accountsList", :minimum, top, 1, false)

    Swagger.validate_param("skip", "accountsList", :minimum, skip, 1, false)

    _ctx = Swagger.Ctx(_api.client, "GET", DataLakeStoreAccountListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.DataLakeStore/accounts", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$skip", skip)  # type Int32
    Swagger.set_param(_ctx.query, "\$select", select)  # type String
    Swagger.set_param(_ctx.query, "\$orderby", orderby)  # type String
    Swagger.set_param(_ctx.query, "\$count", count)  # type Bool
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function accountsList(_api::AccountsApi, subscriptionId::String, api_version::String; filter=nothing, top=nothing, skip=nothing, select=nothing, orderby=nothing, count=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_accountsList(_api, subscriptionId, api_version; filter=filter, top=top, skip=skip, select=select, orderby=orderby, count=count, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function accountsList(_api::AccountsApi, response_stream::Channel, subscriptionId::String, api_version::String; filter=nothing, top=nothing, skip=nothing, select=nothing, orderby=nothing, count=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_accountsList(_api, subscriptionId, api_version; filter=filter, top=top, skip=skip, select=select, orderby=orderby, count=count, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists the Data Lake Store accounts within a specific resource group. The response includes a link to the next page of results, if any.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: filter::String
Param: top::Int32
Param: skip::Int32
Param: select::String
Param: orderby::String
Param: count::Bool
Return: DataLakeStoreAccountListResult
"""
function _swaggerinternal_accountsListByResourceGroup(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, api_version::String; filter=nothing, top=nothing, skip=nothing, select=nothing, orderby=nothing, count=nothing, _mediaType=nothing)
    Swagger.validate_param("top", "accountsListByResourceGroup", :minimum, top, 1, false)

    Swagger.validate_param("skip", "accountsListByResourceGroup", :minimum, skip, 1, false)

    _ctx = Swagger.Ctx(_api.client, "GET", DataLakeStoreAccountListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$skip", skip)  # type Int32
    Swagger.set_param(_ctx.query, "\$select", select)  # type String
    Swagger.set_param(_ctx.query, "\$orderby", orderby)  # type String
    Swagger.set_param(_ctx.query, "\$count", count)  # type Bool
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function accountsListByResourceGroup(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, api_version::String; filter=nothing, top=nothing, skip=nothing, select=nothing, orderby=nothing, count=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_accountsListByResourceGroup(_api, subscriptionId, resourceGroupName, api_version; filter=filter, top=top, skip=skip, select=select, orderby=orderby, count=count, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function accountsListByResourceGroup(_api::AccountsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, api_version::String; filter=nothing, top=nothing, skip=nothing, select=nothing, orderby=nothing, count=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_accountsListByResourceGroup(_api, subscriptionId, resourceGroupName, api_version; filter=filter, top=top, skip=skip, select=select, orderby=orderby, count=count, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates the specified Data Lake Store account information.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: parameters::UpdateDataLakeStoreAccountParameters (required)
Param: api_version::String (required)
Return: DataLakeStoreAccount
"""
function _swaggerinternal_accountsUpdate(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", DataLakeStoreAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function accountsUpdate(_api::AccountsApi, subscriptionId::String, resourceGroupName::String, accountName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsUpdate(_api, subscriptionId, resourceGroupName, accountName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function accountsUpdate(_api::AccountsApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_accountsUpdate(_api, subscriptionId, resourceGroupName, accountName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export accountsCheckNameAvailability, accountsCreate, accountsDelete, accountsEnableKeyVault, accountsGet, accountsList, accountsListByResourceGroup, accountsUpdate


text:
#
# FutureValueModel
#

function initial_state(ctx::CompilerContext{F}) where {F<:FutureValueModel}
    return get_pricing_date(ctx)
end

function lower!(ctx::CompilerContext{F}, rf::RiskFactor, state) :: OptimizingIR.ImmutableValue where {F<:FutureValueModel}
    # in StaticHedging, all risk factors are valued in the pricing date
    assert_at_initial_state(ctx, state)
    return add_input_riskfactor!(ctx, rf)
end

function lower!(ctx::CompilerContext{FutureValueModel{C}}, c::Unit{SpotCurrency{C}}, at::Date) :: OptimizingIR.ImmutableValue where {C<:Currencies.Currency}
    assert_at_initial_state(ctx, at)
    # same currency
    return lower!(ctx, Konst(1.0), at)
end

function lower!(ctx::CompilerContext{FutureValueModel{C1}}, c::Unit{SpotCurrency{C2}}, at::Date) :: OptimizingIR.ImmutableValue where {C1<:Currencies.Currency, C2<:Currencies.Currency}
    assert_at_initial_state(ctx, at)
    return lower!(ctx, c.o, at)
end

function lower!(ctx::CompilerContext{F}, o::ObservableAt, at::Date) :: OptimizingIR.ImmutableValue where {F<:FutureValueModel}
    @assert o.at <= at "Observable date was fixed to ($(o.at)), and can't be valued after date ($(at))"
    return lower!(ctx, o.o, o.at)
end

function lower!(ctx::CompilerContext{F}, c::WhenAt{U}, at::Date) :: OptimizingIR.ImmutableValue where {U<:Unit, F<:FutureValueModel}
    assert_at_initial_state(ctx, at)
    if get_pricing_date(ctx) > c.maturity
        # ignore old cashflows
        return lower!(ctx, Worthless(), at)
    end

    # disable discounted cashflow rule for FutureValueModel
    unit = c.c
    return lower!(ctx, unit, at)
end

# TODO: is this a general rule, or only applies to StaticHedgingModel / FutureValueModel ?
function lower!(ctx::CompilerContext{F}, c::WhenAt{Scale}, state) :: OptimizingIR.ImmutableValue where {F<:FutureValueModel}
    assert_at_initial_state(ctx, state)
    pricing_date = get_pricing_date(ctx)

    mat = c.maturity
    scale_contract = c.c
    scale_observable = scale_contract.s
    scale_inner_contract = scale_contract.c

    if pricing_date > mat
        # ignore old cashflows
        return lower!(ctx, Worthless(), state)
    elseif pricing_date == mat
        # returns the underlying asset value if at maturity
        return lower!(ctx, c.c, state)
    end

    # obtains the observable value at the acquisition date
    arg1 = lower!(ctx, scale_observable, mat) # TODO: this should be a problem if we wanted to apply this as a general rule
    arg2 = lower!(ctx, WhenAt(mat, scale_inner_contract), state)

    # arg1 * arg2
    return add_instruction!(ctx, OP_BIN_SCALAR_MUL, arg1, arg2)
end


text:
using Libdl

let AT_DIR = get(ENV, "AT_DIR", (Sys.iswindows() ? "C:/Program Files/Andor SDK3" :
                                 "/usr/local/andor"))
    hdrname = "atcore.h"
    incdir = ""
    if haskey(ENV, "AT_INCDIR")
        dir = ENV["AT_INCDIR"]
        if !isfile(joinpath(dir, hdrname))
            error("\nDirectory specified by environment variable \"AT_INCDIR\" ",
                  "(\"$dir\") does not contain file Andor SDK header file ",
                  "\"$hdrname\".\n",
                  "Fix the definition of \"AT_INCDIR\" and rebuild.")
        end
        incdir = normpath(dir)
    else
        for dir in (joinpath(AT_DIR, "include"), AT_DIR)
            if isdir(dir) && isfile(joinpath(dir, hdrname))
                incdir = normpath(dir)
                break
            end
        end
        if incdir == ""
            error("\nAndor SDK header file \"$hdrname\" not found.\n",
                  "Define environment variable \"AT_INCDIR\" with the ",
                  "directory containing this file and rebuild.")
        end
    end

    libname = (Sys.iswindows() ? "atcore."*Libdl.dlext :
               "libatcore."*Libdl.dlext)
    libdir = ""
    if haskey(ENV, "AT_LIBDIR")
        dir = normpath(ENV["AT_LIBDIR"])
        if !isfile(joinpath(dir, libname))
            error("\nDirectory specified by environment variable \"AT_LIBDIR\" ",
                  "(\"$dir\") does not contain file Andor SDK library file ",
                  "\"$libname\".\n",
                  "Fix the definition of \"AT_LIBDIR\" and rebuild.")
        end
        libdir = normpath(dir)
    else
        for dir in (joinpath(AT_DIR, "lib"), AT_DIR)
            if isdir(dir) && isfile(joinpath(dir, libname))
                libdir = normpath(dir)
                break
            end
        end
        if libdir == ""
            error("\nAndor SDK library file \"$libname\" not found.\n",
                  "Define environment variable \"AT_LIBDIR\" with the ",
                  "directory containing this file and rebuild.")
        end
    end
    dll = joinpath(libdir, libname)

    # FIXME: Building all require old usb.h to be available.
    #target = (Sys.islinux() ? "all" : "deps.jl")
    target = "deps.jl"

    print(stderr, "\nI will run the following command:\n",
            "    make $target AT_INCDIR=\"$incdir\" ",
            "AT_LIBDIR=\"$libdir\" AT_DLL=\"$dll\"\n\n")
    run(`make $target AT_INCDIR="$incdir" AT_LIBDIR="$libdir" AT_DLL="$dll"`)
end


text:
__precompile__()
module SimpleDirectMediaLayer

    using ColorTypes

    import Base.unsafe_convert
    export  TTF_Init, TTF_OpenFont, TTF_RenderText_Blended, TTF_SizeText

    const depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
    if isfile(depsfile)
        include(depsfile)
    else
        error("SimpleDirectMediaLayer not properly installed. Please run Pkg.build(\"SimpleDirectMediaLayer\") then restart Julia.")
    end

    include("lib/SDL.jl")
    include("lib/SDL_ttf.jl")
    include("lib/SDL_mixer.jl")

    mutable struct SDLWindow
        win::Ptr{Window}
        renderer::Ptr{Renderer}

        function SDLWindow(w,h,title="SDL Window")
            win = CreateWindow(title, Int32(100), Int32(100), Int32(w), Int32(h), UInt32(WINDOW_SHOWN | WINDOW_INPUT_FOCUS))
            SetWindowResizable(win,true)

            renderer = CreateRenderer(win, Int32(-1), UInt32(RENDERER_ACCELERATED | RENDERER_PRESENTVSYNC) )

            new(win, renderer)
        end
    end

    function init()
        GL_SetAttribute(GL_MULTISAMPLEBUFFERS, 4)
        GL_SetAttribute(GL_MULTISAMPLESAMPLES, 4)
        Init(UInt32(INIT_VIDEO))
        TTF_Init()
        Mix_OpenAudio(Int32(22050), UInt16(MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )
    end

    function mouse_position()
        x,y = Int[1], Int[1]
        GetMouseState(pointer(x), pointer(y))
        x[1],y[1]
    end

    function event()
        ev = Event(ntuple(i->UInt8(0),56))
        PollEvent(pointer_from_objref(ev)) == 0 && return nothing

        t = UInt32(0)
        for x in ev._Event[4:-1:1]
            t = t << (sizeof(x)*8)
            t |= x
        end
        evtype = Event(t)
        evtype == nothing && return nothing

        unsafe_load( Ptr{evtype}(pointer_from_objref(ev)) )
    end

end # module


text:
using Plots
using Reproduce
using JLD2
using FileIO

include("plot_utils.jl")

function plot_single_goal_visitation(results; step_size=10)
    p = plot()

    save_file = "plotting/chunlok/generated_plots/single_goal_visitation.svg"

    visits_arr = results[:visit_order]
    episode_length_arr = results[:episode_length]

    visit_perc = get_single_goal_percentage(visits_arr;step_size=step_size)

    for i in 1:4
        plot!(p, 1:step_size:step_size*size(visit_perc)[2], visit_perc[i, :])
    end

    # print(error[3,10])
    savefig(save_file)
    println("Single goal visitation plot saved to $(save_file)")
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using DrWatson
@quickactivate "StatReth"

# %%
"""2H1
Pr(twins|speciesA) = 0.1
Pr(twins|speciesB) = 0.2
Pr(speciesA) = Pr(speciesB) = 0.5
Pr(twins) = Pr(twins|speciesA) * Pr(speciesA) + Pr(twins|speciesB) * Pr(speciesB)
          = 0.1 * 0.5 + 0.2 * 0.5 = 0.15

Now births twins.

Bayes:
Pr(species|twins) = Pr(twins|species) * Pr(species) / Pr(twins)
Pr(panda = speciesA|twins) = 0.1 * 0.5 / 0.15 = 1/3
Pr(panda = speciesB|twins) = 0.2 * 0.5 / 0.15 = 2/3

Pr(twin|panda) = Pr(twins|speciesA) * Pr(panda = speciesA) +
                 Pr(twins|speciesB) * Pr(panda = speciesB)
               = 0.1 * 1/3 + 0.2 * 2/3 = 0.167

"""
# %%
"""2H2
Pr(twins|speciesA) = 0.1
Pr(twins|speciesB) = 0.2
Pr(speciesA) = Pr(speciesB) = 0.5
Pr(twins) = Pr(twins|speciesA) * Pr(speciesA) + Pr(twins|speciesB) * Pr(speciesB)
          = 0.1 * 0.5 + 0.2 * 0.5 = 0.15

Now births twins.

Bayes:
Pr(species|twins) = Pr(twins|species) * Pr(species) / Pr(twins)
Pr(panda = speciesA|twins) = 0.1 * 0.5 / 0.15 = 1/3
"""
# %%
"""2H3
Pr(twins|speciesA) = 0.1  |  Pr(singleinfant|speciesA) = 0.9
Pr(twins|speciesB) = 0.2  |  Pr(singleinfant|speciesB) = 0.8

Now:
Pr(speciesA) = 1/3
Pr(speciesB) = 2/3
Pr(twin|panda) = 0.167

Pr(singleinfant) = Pr(singleinfant|speciesA) * Pr(speciesA) +
                   Pr(singleinfant|speciesB) * Pr(speciesB)
                 = 0.9 * 1/3 + 0.8 * 2/3 = 0.833

Now births single infant.

Bayes:
Pr(species|singleinfant) = Pr(singleinfant|species) * Pr(species) / Pr(singleinfant)
Pr(panda = speciesA|singleinfant) = 0.9 * 1/3 / 0.833 = 0.36
"""
# %%
"""2H4
I understood the test like this:

                          panda is:
                    speciesA    speciesB
test    speciesA      0.8         0.35
says:   speciesB      0.2         0.65

Now the test comes back as speciesA.

We would like to know Pr(panda = speciesA | test = speciesA).

Bayes:
Pr(test = A | panda = A) = 0.8
Pr(panda = A) = Pr(panda = B) = 0.5
Pr(test = A) = 0.8 * 0.5 + 0.35 * 0.5 = 0.575
Pr(test = B) = 0.2 * 0.5 + 0.65 * 0.5 = 0.425

Pr(panda = A | test = A)
            = Pr(test = A | panda = A) * Pr(panda = A) / Pr(test = A)
            = 0.8 * 0.5 / 0.575 = 0.696
Pr(panda = B | test = A)
            = Pr(test = A | panda = B) * Pr(panda = B) / Pr(test = A)
            = 0.35 * 0.5 / 0.425 = 0.304

This will be our new prior:
=> Pr(speciesA) = 0.696
=> Pr(speciesB) = 0.304

Now births twins:
We want to know: Pr(speciesA|twins)

Pr(twins|speciesA) = 0.1  |  Pr(singleinfant|speciesA) = 0.9
Pr(twins|speciesB) = 0.2  |  Pr(singleinfant|speciesB) = 0.8

Bayes:
Pr(speciesA|twins) = Pr(twins|speciesA) * Pr(speciesA) / Pr(twins)
Pr(twins|speciesA) = 0.1
Pr(speciesA) = 0.696
Pr(twins) = 0.1 * 0.696 + 0.2 * 0.304 = 0.130
=> Pr(speciesA|twins) = 0.1 * 0.696 / 0.130 = 0.535
=> Pr(speciesB|twins) = 0.2 * 0.304 / 0.130 = 0.468
(including some rounding errors)

New priors:
=> Pr(speciesA) = 0.535
=> Pr(speciesB) = 0.468

Now births twins:
We want to know: Pr(speciesA|singleinfant)

Bayes:
Pr(speciesA|singleinfant) = Pr(singleinfant|speciesA) * Pr(speciesA) / Pr(singleinfant)
Pr(singleinfant|speciesA) = 0.9
Pr(speciesA) = 0.535
Pr(singleinfant) = 0.9 * 0.535 + 0.8 * 0.468 = 0.856
=> Pr(speciesA|twins) = 0.9 * 0.535 / 0.856 = 0.563
=> Pr(speciesB|twins) = 0.8 * 0.468 / 0.856 = 0.437

"""

# %%
using Statistics
using BenchmarkTools

function genetest(p)
    if p == true  # panda is A
        return rand() < 0.8
    else          # panda is B
        return rand() > 0.65
    end
end

function twinbirth(p)
    if p == true  # panda is A
        return rand() < 0.1
    else          # panda is B
        return rand() < 0.2
    end
end

function singleinfant(p)
    if p == true  # panda is A
        return rand() < 0.9
    else          # panda is B
        return rand() < 0.8
    end
end

f(p) = mean(p[@. genetest(p) & twinbirth(p) & singleinfant(p)])

function g(p)
    countA = count(p) do p
        genetest(p) && twinbirth(p) && singleinfant(p) && p
    end
    countB = count(p) do p
        genetest(p) && twinbirth(p) && singleinfant(p) && !p
    end
    countA / (countA + countB)
end

function h(pandasA)
    countA = 0
    countB = 0
    for p in pandasA
        if genetest(p) && twinbirth(p) && singleinfant(p)
            if p
                countA += 1
            else
                countB += 1
            end
        end
    end
    countA / (countA + countB)
end

n = 10000000
pandasA = rand(Bool, n)

f(pandasA)
g(pandasA)
h(pandasA)


text:
export Scattered1d

abstract type Scattered end

"""
    Scattered1d

A structure for scattering coefficients.

## Fields: Scattered1d
 | **Field** | **Description** |
 |:----------|:----------------|
 | :T        | signal length |
 | :nFilters | number of filters per layer |
 | :Coeff    | scattering coefficients |
"""
struct Scattered1d <: Scattered
    T::Int64 # number of samples in time
    nFilters::Vector{Int64} # number of filters in each layer
    Coeff::Vector{Array{Float64,2}} # scattering coefficients for each layer

    function Scattered1d(F::Vector{FilterBank1d}; subsample::Bool=true)
        nFilters = Vector{Int64}(undef, length(F)+1)
        Coeff = Vector{Array{Float64,2}}(undef, length(F)+1)

        # 0th order scattering coefficients
        nFilters[1] = 1
        T = 2^F[1].N
        Coeff[1] = Array{Float64,2}(undef, T, 1)

        # higher order scattering coefficients
        nVect=1
        for i=1:length(F)
            nFilters[i+1] = length(F[i].Λ)
            if i>=2 && subsample
                # count the number of js in the second layer greater than js in the first layer
                nVect=sum([length([ψ2.j for ψ2 in F[i].Λ if ψ2.j>ψ1.j]) for ψ1 in F[i-1].Λ])
            else
                nVect*=length(F[i].Λ)
            end
            Coeff[i+1] = Array{Float64,2}(undef, T, nVect)
        end

        return new(T, nFilters, Coeff)
    end

    # TODO: Reduce the number of convolutions by counting the js[i]<js[i-1]
    function Scattered1d(F::Vector{FilterBank1dBlock})
        nFilters = Vector{Int64}(undef, length(F)+1)
        Coeff = Vector{Array{Float64,2}}(undef, length(F)+1)

        # 0th order scattering coefficients
        nFilters[1] = 1
        T = 2^F[1].N
        Coeff[1] = Array{Float64,2}(undef, T, 1)

        # higher order scattering coefficients
        nRows=1
        for i=1:length(F)
            nFilters[i+1] = size(F[i].Λ, 2)
            nRows=nRows*size(F[i].Λ, 2)
            Coeff[i+1] = Array{Float64,2}(undef, T, nRows)
        end

        return new(T, nFilters, Coeff)
    end
end


text:
function fwd_prop(X::Array{Float64,1},
                  nodes::Int64,
                  hidden_layer::Int64,
                  z,
                  layer,
                  layer_p1)
    # number of used weights
    j = 0

    # number of covariates
    D = length(X)

    # build first hidden layer
    for node in 1:nodes
        i = j + 1
        j = i + D
        layer[node] = sigmoid(dot(z[i:(j-1)], X) + z[j])
    end

    # build other hidden layers
    if hidden_layer > 1
        for l in 2:hidden_layer
            for node in 1:nodes
                i = j + 1
                j = i + nodes
                layer_p1[node] = sigmoid(dot(z[i:(j-1)], layer) + z[j])
            end

            layer = copy(layer_p1)
        end
    end

    # build output
    i = j + 1
    j = i + nodes
    out = sigmoid(dot(z[i:(j-1)], layer) + z[j])

    return out
end


text:
TA{F, N} = Tracker.TrackedArray{F, N, Array{F, N}}
TR{F} = Tracker.TrackedReal{F}

abstract type Advi end
abstract type VP <: Advi end
abstract type RealSpace <: Advi end
abstract type TranSpace <: Advi end

mutable struct State{T <: Advi, F, A1, A2, A3}
  delta0::A1
  delta1::A1
  sig2::A1
  W::A2
  eta0::A3
  eta1::A3
  v::A1
  H::A2
  alpha::A1 # F won't work, TrackedReals don't work as expected

  State(T::Type, F::Type, A::Type) = new{T, F, A{1}, A{2}, A{3}}()
end

function rsample(s::State{VP, ADVI.MPR{F}, ADVI.MPA{F, 1}, ADVI.MPA{F, 2}, ADVI.MPA{F, 3}};
                 RT::Type=TR{F}, AT::Type=TA{F}) where {F <: AbstractFloat}

  # Or, RT=Float64, AT=Array{Float64}

  real = State(RealSpace, RT, AT)
  tran = State(TranSpace, RT, AT)

  for key in fieldnames(State)
    f = getfield(s, key)
    if typeof(f) <: Array
      # TODO: optimize
      rs = []
      ts = []
      for each_f in f
        println(f)
        r = ADVI.rsample(f)
        t = ADVI.transform(f, r)
        append!(rs, r)
        append!(ts, t)
      end
      setfield!(real, key, rs)
      setfield!(tran, key, ts)

    else
      r = ADVI.rsample(f)
      t = ADVI.transform(f, r)
      setfield!(real, key, r)
      setfield!(tran, key, t)
    end
  end

  return real, tran
end


text:
using Base
using JSON, Unmarshal
using SynchronySDK
using UUIDs

# 1. Get a Synchrony configuration
# Assume that you're running in local directory
cd(joinpath(Pkg.dir("SynchronySDK"),"examples"))
configFile = open("synchronyConfig_Local.json")
configData = JSON.parse(readstring(configFile))
close(configFile)
synchronyConfig = Unmarshal.unmarshal(GraffConfig, configData)

# 2a. User creation and retrieval
# In Synchrony, users are the root of all robots,
# so after creating a user, you can create a bot
# and associate it with data sessions
# Note that we'll create this for you normally, but providing it here to demonstrate
# how a user is associated with an organization.
newUser = UserRequest(synchronyConfig.userId, "NewUser", "email@email.com", "N/A", "Student", "Student", string(uuid4()))
retUser = addUser(synchronyConfig, newUser)
@show retUser
# Now we can get it as well if we want
getUser = getUser(synchronyConfig, newUser.id)
if (JSON.json(retUser) != JSON.json(getUser))
    error("Hmm, users should match")
end

# 2b. Config retrieval
# During runtime, we would skip step 3a and just retrieve the configuration
# for our user. This contains all the parameters required to ingest or retrieve
# data from the system.
runtimeConfig = getUserConfig(synchronyConfig, newUser.id)


text:
using BEASTXMLConstructor

xc = BEASTXMLConstructor

using DataFrames, CSV, Test
### Making sure https://github.com/suchard-group/incomplete_measurements/blob/master/scripts/xml_setup.jl works

function make_mbd(data_path::String, newick_path::String, xml_path::String,
                filename::String; dates_path::String = "")

    df = DataFrame(CSV.File(data_path))

    use_dates = false
    if length(dates_path) > 0
        dates_df = DataFrame(CSV.File(dates_path))
        @assert dates_df[!, :taxon] == df[!, :taxon]
        use_dates = true
    end

    newick = read(newick_path, String)

    taxa = df[!, 1]
    data = Matrix(df[!, 2:end])

    bx = BEASTXMLConstructor.make_residual_xml(data, taxa, newick, chain_length = 100_000)
    if use_dates
        BEASTXMLConstructor.use_dates!(bx)
        xc.set_data_dates(bx, dates_df[!, :date])
    end
    xc.set_screen_logEvery(bx, 100)
    xc.set_file_logEvery(bx, 10)
    xc.set_filename(bx, filename)
    BEASTXMLConstructor.add_MBD_loggables!(bx)

    BEASTXMLConstructor.save_xml(xml_path, bx)
end

n = 10
p = 4
taxa = ["taxon$i" for i = 1:n]
data = randn(n, p)
dates = rand(n)

df = DataFrame()
df.taxon = taxa
for i = 1:p
    df[!, Symbol("trait$i")] = data[:, i]
end

cd(@__DIR__)

data_path = "data.csv"
dates_path = "dates.csv"
newick_path = joinpath("data", "newick.txt")
xml_path = "BEASTXMLConstructor.xml"
dates_xml_path = "xml_dates.xml"
filename = "test"

CSV.write(data_path, df)
CSV.write(dates_path, DataFrame(taxon = taxa, date = dates))

newick = read(newick_path, String)

make_mbd(data_path, newick_path, xml_path, filename)

@test isfile(xml_path)

make_mbd(data_path, newick_path, dates_xml_path, filename, dates_path = dates_path)

@test isfile(dates_xml_path)

### Testing Factor xml
k = 3

# HMC, no shrinkage
bx = BEASTXMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = false)
BEASTXMLConstructor.save_xml("facHMC.xml", bx)
@test isfile("facHMC.xml")

# Gibbs, no shrinkage
bx = BEASTXMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = true)
BEASTXMLConstructor.save_xml("facGibbs.xml", bx)
@test isfile("facGibbs.xml")

# # HMC, shrinkage
# bx = BEASTXMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = true,
#             shrink_loadings = true)
# BEASTXMLConstructor.save_xml("facHMCShrink.xml", bx)
# @test isfile("facHMCShrink.xml")

# # Gibbs, shrinkage
# bx = BEASTXMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = false,
#             shrink_loadings = true)
# BEASTXMLConstructor.save_xml("facGibbsShrink.xml", bx)
# @test isfile("facGibbsShrink.xml")

# # Rotate prior
# fn = "facGibbsRotate.xml"
# bx = BEASTXMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = false,
#             shrink_loadings = true, rotate_prior = true)
# BEASTXMLConstructor.save_xml(fn, bx)
# @test isfile(fn)

################################################################################
## Orthogonal factor analysis
################################################################################

fn = "facOrthogonal.xml"
bx = BEASTXMLConstructor.make_orthogonal_pfa_xml(data, taxa, newick, k)
BEASTXMLConstructor.save_xml(fn, bx)
@test isfile(fn)

################################################################################
## Joint models
################################################################################

n = 10
k = 2
p_res = 3
p_fac = 5
p_diff = k + p_res

taxa = ["taxon_$i" for i = 1:n]

# tree = rtree(taxa)
# newick = writeTopology(tree)

rm = ResidualVarianceModel(p_res)
fm = IntegratedFactorModel(k, p_fac)

jm = JointProcessModel([rm, fm])

dm = DataModel(taxa,
               [randn(n, p_res), randn(n, p_fac)],
               ["trait.res", "trait.fac"])

bx = make_joint_xml(newick, dm, jm)
joint_path = "joint.xml"
# BEASTXMLConstructor.save_xml(joint_path, bx)
# @test isfila(joint_path)
# rm(joint_path)

################################################################################
## Integrate two xml
################################################################################

test_path = joinpath(@__DIR__, "data", "sequence.xml")

bx_seq = BEASTXMLElement(test_path)

taxa = BEASTXMLConstructor.find_element(bx_seq, BEASTXMLConstructor.EmptyDataXMLElement).taxa

k = 2
p = 10
n = length(taxa)
data = randn(n, p)
# newick = writeTopology(rtree(taxa))

bx = BEASTXMLConstructor.make_pfa_xml(data, taxa, newick, k, useHMC = false)
BEASTXMLConstructor.save_xml("facGibbs.xml", bx)
@test isfile("facGibbs.xml")

BEASTXMLConstructor.merge_xml!(bx, bx_seq)
# xml = BEASTXMLConstructor.make_xml(bx)
BEASTXMLConstructor.save_xml("merge.xml", bx);
# print(xml)

################################################################################
## Sampled factor model with integrated factors sampler
################################################################################
k = 2
p = 10
n = length(taxa)
data = randn(n, p)
# newick = writeTopology(rtree(taxa))

bx = BEASTXMLConstructor.make_sampled_pfa_xml(data, taxa, newick, k)
save_xml("facSampledHMC.xml", bx, change_filename=true)

################################################################################
## Latent liability
################################################################################

states = ones(Int, p)

states[p - 1] = 2
data[:, p - 1] .= rand(0:1, n)

states[p] = 4
data[:, p] .= rand(0:3, n)

bx = make_pfa_xml(data, taxa, newick, k)
BEASTXMLConstructor.add_latent_liability(bx, ["traits"], states)
save_xml("liability.xml", bx)


text:
# ODE SIR/SICR models
# variables are concentrations, X/N
# Z = 1-S

function sir(T,N,beta,sigmaR,sigmaD,a,b)

	dI = zeros(Int,T)
	dR = zeros(Int,T)
	dD = zeros(Int,T)
	I = zeros(Int,T)

	S = N
	I[1] = 10
	for i = 2:T
		d = Binomial(S,beta*S*I[i-1]/N^2)
		dI[i] = rand(d)
		d = Binomial(I[i-1],sigmaR*I[i-1]/N)
		dR[i] = rand(d)
		d=Gamma(a,b)
		ind = round(Int,rand(d))
		ind = max(i-1-ind,1)
		d = Binomial(I[i-1],sigmaD*I[ind]/N)
		dD[i] = rand(d)
		I[i] = I[i-1] + dI[i] - dR[i] - dD[i]
		S -=  dI[i]
	end

	return dI, dR, dD

end

function sicr(T,N,beta,sigmaI,sigmaC,sigmaR,sigmaD,a,b)

	dI = zeros(Int,T)
	dC = zeros(Int,T)
	dR = zeros(Int,T)
	dD = zeros(Int,T)
	dC = zeros(Int,T)
	C = zeros(Int,T)
	S = N
	I = 10
	for i = 2:T
		d = Binomial(S,beta*S*I/N^2)
		dIp = rand(d)
		d = Binomial(I,sigmaI*I/N)
		dIm = rand(d)
		d = Binomial(I,sigmaC*I/N)
		dC[i] = rand(d)
		d = Binomial(C[i-1],sigmaR*C[i-1]/N)
		dR[i] = rand(d)
		d=Gamma(a,b)
		ind = round(Int,rand(d))
		ind = max(i-1-ind,1)
		d = Binomial(C[i-1],sigmaD*C[ind]/N)
		dD[i] = rand(d)
		C[i] = C[i-1] + dC[i] - dR[i] - dD[i]
		I += dIp - dIm - dC[i]
		S -=  dIp
	end

	return dC, dR, dD

end

function sicr(T,N,beta,sigmaI,sigmaC,sigmaR,sigmaD)

	dI = zeros(Int,T)
	dC = zeros(Int,T)
	dR = zeros(Int,T)
	dD = zeros(Int,T)
	dC = zeros(Int,T)
	C =1
	R = 0
	S = N
	I = 10
	for i = 2:T
		d = Binomial(S,beta*S*I/N^2)
		dIp = rand(d)
		d = Binomial(I,sigmaI*I/N)
		dIm = rand(d)
		d = Binomial(I,sigmaC*I/N)
		dC[i] = rand(d)
		if C > 0
			d = Binomial(C,sigmaR*C/N)
			dR[i] = rand(d)
		else
			dC[i] = 0
		end
		if R > 0
			d = Binomial(R,sigmaD*R/N)
			dD[i] = rand(d)
		else
			dD[i] =0
		end
		C +=  dC[i] - dR[i] - dD[i]
		R += dR[i] - dD[i]
		I += dIp - dIm - dC[i]
		S -=  dIp
	end

	return dC, dR, dD

end

max_stars_repo_name: nih-niddk-mbs/covid-sicr

max_issues_repo_name: nih-niddk-mbs/covid-sicr

max_forks_repo_name: nih-niddk-mbs/covid-sicr


text:
using Pkg

function main()
    if length(ARGS) != 1
        println("Usage: $PROGRAM_FILE <package name>")
        exit()
    end
    tgt = ARGS[1]
    println("Checking packages which requires $(tgt)...")
    pkgs = search(tgt)
    println("$(length(pkgs)) packages found.")
    for ipkg in pkgs
        println(" - $(ipkg)")
    end
end

function search(tgt)
    cont = Array{String, 1}()
    deps = Pkg.dependencies()
    for (key, value) in deps
        name  = value.name
        ideps = value.dependencies
        tgt in keys(ideps) && push!(cont, name)
    end
    return cont
end

main()

max_stars_repo_name: yoshphys/julia-utility

max_issues_repo_name: yoshphys/julia-utility

max_forks_repo_name: yoshphys/julia-utility