text:
function textDocument_hover_request(params::TextDocumentPositionParams, server::LanguageServerInstance, conn)
    doc = getdocument(server, URI2(params.textDocument.uri))
    x = get_expr1(getcst(doc), get_offset(doc, params.position))
    x isa EXPR && CSTParser.isoperator(x) && resolve_op_ref(x, server)
    documentation = get_hover(x, "", server)
    documentation = get_closer_hover(x, documentation)
    documentation = get_fcall_position(x, documentation)
    documentation = sanitize_docstring(documentation)

    return isempty(documentation) ? nothing : Hover(MarkupContent(documentation), missing)
end

function get_hover(x, documentation::String, server) documentation end

function get_hover(x::EXPR, documentation::String, server)
    if (CSTParser.isidentifier(x) || CSTParser.isoperator(x)) && StaticLint.hasref(x)
        r = refof(x)
        documentation = if r isa StaticLint.Binding
            get_hover(r, documentation, server)
        elseif r isa SymbolServer.SymStore
            get_hover(r, documentation, server)
        else
            documentation
        end
    end
    return documentation
end

function get_hover(b::StaticLint.Binding, documentation::String, server)
    if b.val isa StaticLint.Binding
        documentation = get_hover(b.val, documentation, server)
    elseif b.val isa EXPR
        if CSTParser.defines_function(b.val) || CSTParser.defines_datatype(b.val)
            documentation = get_func_hover(b, documentation, server)
            for r in b.refs
                method = StaticLint.get_method(r)
                if method isa EXPR
                    documentation = get_preceding_docs(method, documentation)
                    if CSTParser.defines_function(method)
                        documentation = string(ensure_ends_with(documentation), "```julia\n", Expr(CSTParser.get_sig(method)), "\n```\n")
                    elseif CSTParser.defines_datatype(method)
                        documentation = string(ensure_ends_with(documentation), "```julia\n", Expr(method), "\n```\n")
                    end
                elseif method isa SymbolServer.SymStore
                    documentation = get_hover(method, documentation, server)
                end
            end
        else
            documentation = try
                documentation = if binding_has_preceding_docs(b)
                    string(documentation, Expr(parentof(b.val).args[3]))
                elseif const_binding_has_preceding_docs(b)
                    string(documentation, Expr(parentof(parentof(b.val)).args[3]))
                else
                    documentation
                end
                documentation = string(documentation, "```julia\n", prettify_expr(Expr(b.val)), "\n```\n")
            catch err
                @error "get_hover failed to convert Expr" exception = (err, catch_backtrace())
                throw(LSHoverError(string("get_hover failed to convert Expr")))
            end
        end
    elseif b.val isa SymbolServer.SymStore
        documentation = get_hover(b.val, documentation, server)
    end
    return documentation
end

function prettify_expr(ex::Expr)
    if ex.head === :kw && length(ex.args) == 2
        string(ex.args[1], " = ", ex.args[2])
    else
        string(ex)
    end
end

prettify_expr(ex) = string(ex)

# print(io, x::SymStore) methods are defined in SymbolServer
function get_hover(b::SymbolServer.SymStore, documentation::String, server)
    if !isempty(b.doc)
        documentation = string(documentation, b.doc, "\n")
    end
    documentation = string(documentation, "```julia\n", b, "\n```")
end

function get_hover(f::SymbolServer.FunctionStore, documentation::String, server)
    if !isempty(f.doc)
        documentation = string(documentation, f.doc, "\n")
    end

    documentation = string(documentation, "`$(f.name)` is a `Function`.\n")
    nm = length(f.methods)
    documentation = string(documentation, "**$(nm)** method", nm == 1 ? "" : "s", " for function ", '`', f.name, '`', '\n')
    for m in f.methods
        io = IOBuffer()
        print(io, m.name, "(")
        nsig = length(m.sig)
        for (i, sig) = enumerate(m.sig)
            if sig[1] ≠ Symbol("#unused#")
                print(io, sig[1])
            end
            print(io, "::", sig[2])
            i ≠ nsig && print(io, ", ")
        end
        print(io, ")")
        sig = String(take!(io))

        text = replace(string(m.file, ':', m.line), "\\" => "\\\\")
        link = text

        if server.clientInfo !== missing
            if occursin("code", lowercase(server.clientInfo.name)) && isabspath(m.file)
                link = string(filepath2uri(m.file), "#", m.line)
            end
        end

        documentation = string(documentation, "- `$(sig)` in `$(m.mod)` at [$(text)]($(link))", '\n')
    end
    return documentation
end

get_func_hover(x, documentation, server) = documentation
get_func_hover(x::SymbolServer.SymStore, documentation, server) = get_hover(x, documentation, server)

function get_preceding_docs(expr::EXPR, documentation)
    if expr_has_preceding_docs(expr)
        string(documentation, Expr(parentof(expr).args[3]))
    elseif is_const_expr(parentof(expr)) && expr_has_preceding_docs(parentof(expr))
        string(documentation, Expr(parentof(parentof(expr)).args[3]))
    else
        documentation
    end
end

ensure_ends_with(s, c = "\n") = endswith(s, c) ? s : string(s, c)

binding_has_preceding_docs(b::StaticLint.Binding) = expr_has_preceding_docs(b.val)

function const_binding_has_preceding_docs(b::StaticLint.Binding)
    p = parentof(b.val)
    is_const_expr(p) && expr_has_preceding_docs(p)
end

expr_has_preceding_docs(x) = false
expr_has_preceding_docs(x::EXPR) = is_doc_expr(parentof(x))

is_const_expr(x) = false
is_const_expr(x::EXPR) = headof(x) === :const

is_doc_expr(x) = false
function is_doc_expr(x::EXPR)
    return CSTParser.ismacrocall(x) &&
        length(x.args) == 4 &&
        headof(x.args[1]) === :globalrefdoc &&
        CSTParser.isstring(x.args[3])
end

get_fcall_position(x, documentation, visited=nothing) = documentation

function get_fcall_position(x::EXPR, documentation, visited=EXPR[])
    if x in visited                                      # TODO: remove
        throw(LSInfiniteLoop("Possible infinite loop.")) # TODO: remove
    else                                                 # TODO: remove
        push!(visited, x)                                # TODO: remove
    end                                                  # TODO: remove
    if parentof(x) isa EXPR
        if CSTParser.iscall(parentof(x))
            call_counts = StaticLint.call_nargs(parentof(x))
            call_counts[1] < 5 && return documentation
            arg_i = 0
            for (i, arg) = enumerate(parentof(x))
                if arg == x
                    arg_i = div(i - 1, 2)
                    break
                end
            end
            arg_i == 0 && return documentation
            fname = CSTParser.get_name(parentof(x))
            if StaticLint.hasref(fname) &&
                (refof(fname) isa StaticLint.Binding && refof(fname).val isa EXPR && CSTParser.defines_struct(refof(fname).val) && StaticLint.struct_nargs(refof(fname).val)[1] == call_counts[1])
                dt_ex = refof(fname).val
                args = dt_ex.args[3]
                args.args === nothing || arg_i > length(args.args) && return documentation
                _fieldname = CSTParser.str_value(CSTParser.get_arg_name(args.args[arg_i]))
                documentation = string("Datatype field `$_fieldname` of $(CSTParser.str_value(CSTParser.get_name(dt_ex)))", "\n", documentation)
            elseif StaticLint.hasref(fname) && (refof(fname) isa SymbolServer.DataTypeStore || refof(fname) isa StaticLint.Binding && refof(fname).val isa SymbolServer.DataTypeStore)
                dts = refof(fname) isa StaticLint.Binding ? refof(fname).val : refof(fname)
                if length(dts.fieldnames) == call_counts[1] && arg_i <= length(dts.fieldnames)
                    documentation = string("Datatype field `$(dts.fieldnames[arg_i])`", "\n", documentation)
                end
            else
                documentation = string("Argument $arg_i of $(call_counts[1]) in call to `", CSTParser.str_value(fname), "`\n", documentation)
            end
            return documentation
        else
            return get_fcall_position(parentof(x), documentation, visited)
        end
    end
    return documentation
end

# info on what expression the current token (e.g. a ], ), `end`, etc.)
get_closer_hover(x, documentation) = documentation
function get_closer_hover(x::EXPR, documentation)
    if parentof(x) isa EXPR
        if headof(x) === :END
            if headof(parentof(x)) === :function
                documentation = string(documentation, "Closes function definition for `", Expr(CSTParser.get_sig(parentof(x))), "`\n")
            elseif CSTParser.defines_module(parentof(x)) && length(parentof(x).args) > 1
                documentation = string(documentation, "Closes module definition for `", Expr(parentof(x).args[2]), "`\n")
            elseif CSTParser.defines_struct(parentof(x))
                documentation = string(documentation, "Closes struct definition for `", Expr(CSTParser.get_sig(parentof(x))), "`\n")
            elseif headof(parentof(x)) === :for && length(parentof(x).args) > 2
                documentation = string(documentation, "Closes for-loop expression over `", Expr(parentof(x).args[2]), "`\n")
            elseif headof(parentof(x)) === :while && length(parentof(x).args) > 2
                documentation = string(documentation, "Closes while-loop expression over `", Expr(parentof(x).args[2]), "`\n")
            else
                documentation = "Closes `$(headof(parentof(x)))` expression."
            end
        elseif headof(x) === :RPAREN
            if CSTParser.iscall(parentof(x)) && length(parentof(x).args) > 0
                documentation = string(documentation, "Closes call of ", Expr(parentof(x).args[1]), "\n")
            end
        end
    end
    return documentation
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

import FMIImport.FMICore: fmi3VariableNamingConventionFlat

myFMU = fmi3Load("BouncingBall", "ModelicaReferenceFMUs", "0.0.14")

@test fmi3GetVersion(myFMU) == "3.0-beta.5"

@test fmi3GetModelName(myFMU) == "BouncingBall"
@test fmi3GetVariableNamingConvention(myFMU) == fmi3VariableNamingConventionFlat
@test fmi3IsCoSimulation(myFMU) == true
@test fmi3IsModelExchange(myFMU) == true
# TODO scheduledExecution
@test fmi3GetInstantiationToken(myFMU) == "{8c4e810f-3df3-4a00-8276-176fa3c9f003}" # TODO update
@test fmi3GetGenerationTool(myFMU) == "Reference FMUs (0278b1e)"
@test fmi3GetGenerationDateAndTime(myFMU) == "[Unknown generation date and time]"
@test fmi3GetNumberOfEventIndicators(myFMU) == 1
@test fmi3CanGetSetState(myFMU) == true
@test fmi3CanSerializeFMUState(myFMU) == true
@test fmi3ProvidesDirectionalDerivatives(myFMU) == false
@test fmi3ProvidesAdjointDerivatives(myFMU) == false

@test fmi3GetDefaultStartTime(myFMU.modelDescription) ≈ 0.0
@test fmi3GetDefaultStopTime(myFMU.modelDescription) ≈ 3.0
@test fmi3GetDefaultTolerance(myFMU.modelDescription) ≈ 1e-4
@test fmi3GetDefaultStepSize(myFMU.modelDescription) === 0.01

fmi3Unload(myFMU)

ext: jl


text:
"""
Computes r̂ across a set of chains from different samplers.
* `schains`: a vector of chains from different samplers
* `csr̂`: A DataFrame containing cross sampler r̂. This is concatonated to the
results DataFrame
"""
function cross_samplerRhat!(schains,csr̂;kwargs...)
    schains = Chains.(schains)
    schains = standardizeNames.(schains)
    chains = reduce(chainscat,schains)
    chains = removeBurnin(chains;kwargs...)
    df = MCMCChains.describe(chains)[1].df
    sort!(df)
    parms = sort!(chains.name_map.parameters)
    values = df[!,:r_hat]
    N = length(schains)
    newDF = DataFrame()
    for (p,v) in zip(parms,values)
        colname = createName(p,:sampler_rhat)
        V = fill(v,N)
        newDF[!,colname] = V
    end
    return vcat(csr̂,newDF,cols=:union)
end

"""
Adds keyword arguments to the results DataFrame
* `df`: DataFrame containing benchmark results for single iteration
* `kwargs`: keyword arguments
"""
function addKW!(df;kwargs...)
    for (k,v) in pairs(kwargs)
        df[!,k] = [v]
    end
end

"""
Adds performance metrics to benchmark results, which include runtime,
memory allocations in MB, garbage collection time, percent of time spent in garbage collection,
and the number of memory allocations
* `df`: the dataframe to which performance metrics are added
* `p`: a collection of performance metrics including run time,
memory allocations and garbage collection time
"""
function addPerformance!(df,p)
    df[!,:time] = [p[2]]
    df[!,:megabytes] = [p[3]/1e6]
    df[!,:gctime] = [p[4]]
    df[!,:gcpercent] = [p[4]/p[2]]
    df[!,:allocations] = [p[5].poolalloc + p[5].malloc]
end

function toDict(data)
    return Dict(string(k)=>v for (k,v) in pairs(data))
end

function removeBurnin(chn;Nadapt,kwargs...)
    return chn[(Nadapt+1):end,:,:]
end

function savechain!(s,chains,performance;savechain=false,kwargs...)
    if savechain
        k = s.name
        push!(chains[k],performance[1])
    end
end

"""
Adds chain summary (e.g. rhat,ess) to newDF for each parameter.
* `newDF`: dataframe that collects results on an iteration
* `chn`: chain for given iteration
* `df`: df of chain results
* `col`: name of column

e.g. If col = :ess, and parameters are mu and sigma, the new columns
will be mu_ess and sigma_ess and will contain their respective ess values
"""
function addChainSummary!(newDF,chn,df,col)
    parms = sort!(chn.name_map.parameters)
    values = df[!,col]
    for (p,v) in zip(parms,values)
        colname = createName(p,col)
        newDF[!,colname] = [v]
    end
end

"""
Effective Sample Size per second
* `newDF`: new DataFrame to be added to results DataFrame
* `chain`: MCMC chain
* `df`: DataFrame containing posterior summaries
* `performance`: a collection of performance metrics including run time,
memory allocations and garbage collection time
"""
function addESStime!(newDF,chn,df,performance)
    parms = sort!(chn.name_map.parameters)
    values = df[!,:ess]/performance[2]
    for (p,v) in zip(parms,values)
        colname = createName(p,"ess_ps")
        newDF[!,colname] = [v]
    end
end

"""
Add highest probability density interval for each parameter
* `newDF`: new DataFrame to be added to results DataFrame
* `chain`: MCMC chain
"""
function addHPD!(newDF,chain)
    h = hpd(chain)
    for r in eachrow(h.df)
        p,lb,ub = r
        colname = createName(string(p),"hdp_lb")
        newDF[!,colname] = [lb]
        colname = createName(string(p),"hdp_ub")
        newDF[!,colname] = [ub]
    end
end

"""
Add mean for each parameter
* `newDF`: new DataFrame to be added to results DataFrame
* `df`: DataFrame containing posterior summaries
"""
function addMeans!(newDF,df)
    for r in eachrow(df)
        p=r[:parameters]
        v = r[:mean]
        colname = createName(string(p),"mean")
        newDF[!,colname] = [v]
    end
end

function standardizeNames(chain)
    nms = names(chain)
    newNames = standardizeName.(nms)
    newChain = Chains(chain.value.data,newNames)
    return newChain
end

function standardizeName(p)
    if occursin(".",p)
        s = split(p,".")
        p = string(s[1],"[",s[2],"]")
    end
    return p
end

function createName(p,col)
    p = standardizeName(p)
    return Symbol(string(p,"_",col))
end

function gettype(s)
    return Symbol(typeof(s).name)
end

initChains(s::T) where {T<:MCMCSampler} = initChains((s,))

"""
Initializes NamedTuple of chains
* s: sampler object
"""
function initChains(s)
    names = map(x->x.name,s)
    v = map(x->Chains[],1:length(s))
    return NamedTuple{names}(v)
end

function Chains(chain::Chains)
    parms = String.(chain.name_map.parameters)
    sort!(parms)
    v=chain[:,parms,:].value.data
    return Chains(v,parms)
end

"""
Run all benchmarks
"""
function run_all_benchmarks()
    path = pathof(MCMCBenchmarks)
    models = [
        "Gaussian",
        "Hierarchial_Poisson",
        "LBA",
        "Linear_Regression",
        "SDT",
        "Autodiff",
    ]
    res = map(models) do m
        @eval module $(Symbol("Test_", m))
            include(joinpath($path, "../../Examples/" * $m * "/" * $m * "_Example.jl"))
        end
    end
 end

 function run_all_benchmarks(models)
     res = map(models) do m
         @eval module $(Symbol("Test_", m))
             include($m)
         end
     end
  end


text:
export induce, InducedMap

struct InducedMap{M<:AbstractIntervalMap, D<:Domain, R<:Domain, HD, I} <: AbstractMarkovMap{D,R}
    he::HofbauerExtension{I,HD,M}
    domain::D
    rangedomain::R
    d_ind::Int
    r_ind::Int

    function InducedMap(he::HofbauerExtension{I,HD,M}, d::D, rd::R, d_ind::Int, r_ind::Int) where
            {M<:AbstractIntervalMap,D<:Domain,R<:Domain,HD,I}
    @assert d == convert(Domain,hdomains(he)[d_ind])
    @assert rd == convert(Domain,hdomains(he)[r_ind])
    new{M,D,R,HD,I}(he,d,rd,d_ind,r_ind)
    end
end
function InducedMap(he::HofbauerExtension, domain::Domain, rangedomain::Domain)
    @assert domain ∈ domains.(he)
    @assert rangedomain ∈ domains.(he)
    InducedMap(he, domain, rangedomain, something(findfirst(isequal(domain),domains.(he))),
        something(findfirst(isequal(rangedomain),domains.(he))))
end

InducedMap(he::HofbauerExtension,d_ind::Integer,r_ind::Integer) =
    InducedMap(he,convert(Domain,hdomains(he)[d_ind]),convert(Domain,hdomains(he)[r_ind]),d_ind,r_ind)
InducedMap(he::HofbauerExtension,d,r) = InducedMap(he,convert(Domain,d),convert(Domain,r))
InducedMap(he::HofbauerExtension,d) = InducedMap(he,d,d)
InducedMap(he::HofbauerExtension) = InducedMap(he,1,1)

InducedMap(m::AbstractIntervalMap;maxdepth=100,forcereturn=true) = InducedMap(hofbauerextension(m;maxdepth=maxdepth,forcereturn=forcereturn))
InducedMap(m::AbstractIntervalMap,d;maxdepth=100,forcereturn=trues(d)) = InducedMap(hofbauerextension(m,d;maxdepth=maxdepth,forcereturn=forcereturn),d)
#length(d::Interval) not defined so:
InducedMap(m::AbstractIntervalMap,d::Union{Domain,ClosedInterval};maxdepth=100,forcereturn=true) = InducedMap(hofbauerextension(m,convert(Domain,d);maxdepth=maxdepth,forcereturn=forcereturn),d)

induce(m::AbstractIntervalMap,args...) = InducedMap(m, args...)

Base.show(io::IO, i::InducedMap) = print(io, "Induced map on $(i.domain)→$(i.rangedomain) of $(i.he.m)")

function (im::InducedMap)(x;maxiter=10000)
    @assert x ∈ domain(im)
    m = getmap(im.he)
    graphind = im.d_ind
    for i = 1:maxiter
        bind = getbranchind(m,x)
        x = maphb(branches(m)[bind],x)
        graphinds = searchsorted(im.he.fedgelist[graphind],HofbauerEdge(0,0,bind,false),by=branchind)
        isempty(graphinds) && error("reached end of generated Hofbauer extension at depth $(hdomains(im.he)[graphind].depth)")
        for g in graphinds
            dst_ind = dst(im.he.fedgelist[graphind][g])
            if x ∈ convert(Domain,hdomains(im.he)[dst_ind])
                graphind = dst_ind
                break
            end
        end
        graphind == im.r_ind && return x
    end
    error("reached max number of iterations: $maxiter")
end

function mapP(im::InducedMap,x;maxiter=10000)
    @assert x ∈ domain(im)
    m = getmap(im.he)
    dx = one(typeof(x))
    p = HofbauerPoint(x,im.d_ind)
    for i = 1:maxiter
        (p,ds) = mapP(im.he,p)
        dx *= ds
        p.graphind == im.r_ind && return (x, dx)
    end
    error("reached max number of iterations: $maxiter")
end
mapD(im::InducedMap,x;maxiter=10000) = mapP(im,x;maxiter=maxiter)[2]

# transferfunction
struct TransferBSFunction{F}
    f::F
end
(q::TransferBSFunction)(v,dvdx,n) = dvdx*q.f(v)
transferfunction(x,im::InducedMap,f) = BackSum(im)(TransferBSFunction(f),true)(x)

struct TransferIntBSFunction{F}
    f::F
end
(q::TransferIntBSFunction)(v,dvdx,n) = q.f(v)
function transferfunction_int(x,y,im::InducedMap,f)
    SQ = BackSum(im)(TransferIntBSFunction(cumsum(f)),true)
    SQ(y) - SQ(x)
end


text:
# Autogenerated wrapper script for OpenEXR_jll for x86_64-apple-darwin
export libIex, libIlmThread, libOpenEXR, libOpenEXRCore, libOpenEXRUtil

using Imath_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("OpenEXR")
JLLWrappers.@declare_library_product(libIex, "@rpath/libIex-3_1.30.dylib")
JLLWrappers.@declare_library_product(libIlmThread, "@rpath/libIlmThread-3_1.30.dylib")
JLLWrappers.@declare_library_product(libOpenEXR, "@rpath/libOpenEXR-3_1.30.dylib")
JLLWrappers.@declare_library_product(libOpenEXRCore, "@rpath/libOpenEXRCore-3_1.30.dylib")
JLLWrappers.@declare_library_product(libOpenEXRUtil, "@rpath/libOpenEXRUtil-3_1.30.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Imath_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libIex,
        "lib/libIex-3_1.30.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libIlmThread,
        "lib/libIlmThread-3_1.30.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOpenEXR,
        "lib/libOpenEXR-3_1.30.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOpenEXRCore,
        "lib/libOpenEXRCore-3_1.30.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOpenEXRUtil,
        "lib/libOpenEXRUtil-3_1.30.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DataFrames
using CSV
using HTTP
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFramesMeta
using GLM
function create_grids()

    function xgrid(theta,xval)
        N      = length(xval)
        xub    = vcat(xval[2:N],Inf)
        xtran1  = zeros(N,N)
        xtran1c = zeros(N,N)
        lcdf   = zeros(N)
        for i=1:length(xval)
            xtran1[:,i]   = (xub[i] .>= xval).*(1 .- exp.(-theta*(xub[i] .- xval)) .- lcdf)
            lcdf        .+= xtran1[:,i]
            xtran1c[:,i] .+= lcdf
        end
        return xtran1, xtran1c
    end

    zval   = collect([.25:.01:1.25]...)
    zbin   = length(zval)
    xval   = collect([0:.125:25]...)
    xbin   = length(xval)
    tbin   = xbin*zbin
    xtran  = zeros(tbin,xbin)
    xtranc = zeros(xbin,xbin,xbin)
    for z=1:zbin
        xtran[(z-1)*xbin+1:z*xbin,:],xtranc[:,:,z] = xgrid(zval[z],xval)
    end

    return zval,zbin,xval,xbin,xtran
end
#Question 1

    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
    df = CSV.read(HTTP.get(url).body,DataFrame)
    # create bus id variable
df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#question 2

theta=glm(@formula(Y ~ Odometer * Odometer^2 * RouteUsage * RouteUsage^2 * Branded * time * time^2), df_long, Binomial(), LogitLink())

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Q3:generate state transition matrices
#:::::::::::::::::::::::::::::::::::::::::::::::::::
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body,DataFrame)
Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
X = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
Z = Vector(df[:,:RouteUsage])
B = Vector(df[:,:Branded])
N = size(Y,1)
T = size(Y,2)
Xstate = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
Zstate = Vector(df[:,:Zst])
zval,zbin,xval,xbin,xtran = create_grids()

#create data frame#
log_DF=DataFrame()
log_DF.Odometer=kron(ones(zbin),xval)
log_DF.RouteUsage =kron(ones(xbin),zval)
log_DF.Branded=zeros(size())
log_DF.time=zeros(size())
#function
function values(Xstate,Zstate,xtran,xbin,zbin,xval)
    FV=zeros(size(xtran,1),2,T+1)
   #Nested loops
    for t=2:T
        for b=0:1
FVT1[i,t] = (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,B[i]+1,t+1]
df_long = @transform(df_long,fv = fvt1)
theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded),
df_long, Binomial(), LogitLink(), offset=df_long.fv)


text:
# Autogenerated wrapper script for Clipper_jll for i686-w64-mingw32-cxx03
export libcclipper

JLLWrappers.@generate_wrapper_header("Clipper")
JLLWrappers.@declare_library_product(libcclipper, "libcclipper.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcclipper,
        "bin\\libcclipper.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Calculate the vorticity on a single node using the derivative matrix from the polynomial basis of
# a DGSEM solver. `u` is the solution on the whole domain.
# This function is used for calculating acoustic source terms for coupled Euler-acoustics
# simulations.
function calc_vorticity_node(u, mesh::TreeMesh{2}, equations::CompressibleEulerEquations2D,
                             dg::DGSEM, cache, i, j, element)
  @unpack derivative_matrix = dg.basis

  v2_x = zero(eltype(u)) # derivative of v2 in x direction
  for ii in eachnode(dg)
    rho, _, rho_v2 = get_node_vars(u, equations, dg, ii, j, element)
    v2 = rho_v2 / rho
    v2_x = v2_x + derivative_matrix[i, ii] * v2
  end

  v1_y = zero(eltype(u)) # derivative of v1 in y direction
  for jj in eachnode(dg)
    rho, rho_v1 = get_node_vars(u, equations, dg, i, jj, element)
    v1 = rho_v1 / rho
    v1_y = v1_y + derivative_matrix[j, jj] * v1
  end

  return (v2_x - v1_y) * cache.elements.inverse_jacobian[element]
end

# Convenience function for calculating the vorticity on the whole domain and storing it in a
# preallocated array
function calc_vorticity!(vorticity, u, mesh::TreeMesh{2}, equations::CompressibleEulerEquations2D,
                         dg::DGSEM, cache)
  @threaded for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      vorticity[i, j, element] = calc_vorticity_node(u, mesh, equations, dg, cache, i, j, element)
    end
  end

  return nothing
end

end # muladd

# From here on, this file contains specializations of DG methods on the
# TreeMesh2D to the compressible Euler equations.
#
# The specialized methods contain relatively verbose and ugly code in the sense
# that we do not use the common "pointwise physics" interface of Trixi.jl.
# However, this is currently (November 2021) necessary to get a significant
# speed-up by using SIMD instructions via LoopVectorization.jl.
#
# TODO: SIMD. We could get even more speed-up if we did not need to permute
#             array dimensions, i.e., if we changed the basic memory layout.
#
# We do not wrap this code in `@muladd begin ... end` block. Optimizations like
# this are handled automatically by LoopVectorization.jl.

# We specialize on `PtrArray` since these will be returned by `Trixi.wrap_array`
# if LoopVectorization.jl can handle the array types. This ensures that `@turbo`
# works efficiently here.
@inline function split_form_kernel!(_du::PtrArray, u_cons::PtrArray,
                                    element, mesh::TreeMesh{2},
                                    nonconservative_terms::Val{false},
                                    equations::CompressibleEulerEquations2D,
                                    volume_flux::typeof(flux_shima_etal_turbo),
                                    dg::DGSEM, cache, alpha)
  @unpack derivative_split = dg.basis

  # Create a temporary array that will be used to store the RHS with permuted
  # indices `[i, j, v]` to allow using SIMD instructions.
  # `StrideArray`s with purely static dimensions do not allocate on the heap.
  du = StrideArray{eltype(u_cons)}(undef,
    (ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))...,
     StaticInt(nvariables(equations))))

  # Convert conserved to primitive variables on the given `element`.
  u_prim = StrideArray{eltype(u_cons)}(undef,
    (ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))...,
     StaticInt(nvariables(equations))))

  @turbo for j in eachnode(dg), i in eachnode(dg)
    rho    = u_cons[1, i, j, element]
    rho_v1 = u_cons[2, i, j, element]
    rho_v2 = u_cons[3, i, j, element]
    rho_e  = u_cons[4, i, j, element]

    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))

    u_prim[i, j, 1] = rho
    u_prim[i, j, 2] = v1
    u_prim[i, j, 3] = v2
    u_prim[i, j, 4] = p
  end

  # x direction
  # At first, we create new temporary arrays with permuted memory layout to
  # allow using SIMD instructions along the first dimension (which is contiguous
  # in memory).
  du_permuted = StrideArray{eltype(u_cons)}(undef,
    (StaticInt(nnodes(dg)), StaticInt(nnodes(dg)),
     StaticInt(nvariables(equations))))

  u_prim_permuted = StrideArray{eltype(u_cons)}(undef,
    (StaticInt(nnodes(dg)), StaticInt(nnodes(dg)),
     StaticInt(nvariables(equations))))

  @turbo for v in eachvariable(equations),
             j in eachnode(dg),
             i in eachnode(dg)
    u_prim_permuted[j, i, v] = u_prim[i, j, v]
  end
  fill!(du_permuted, zero(eltype(du_permuted)))

  # Next, we basically inline the volume flux. To allow SIMD vectorization and
  # still use the symmetry of the volume flux and the derivative matrix, we
  # loop over the triangular part in an outer loop and use a plain inner loop.
  for i in eachnode(dg), ii in (i+1):nnodes(dg)
    @turbo for j in eachnode(dg)
      rho_ll = u_prim_permuted[j, i, 1]
      v1_ll  = u_prim_permuted[j, i, 2]
      v2_ll  = u_prim_permuted[j, i, 3]
      p_ll   = u_prim_permuted[j, i, 4]

      rho_rr = u_prim_permuted[j, ii, 1]
      v1_rr  = u_prim_permuted[j, ii, 2]
      v2_rr  = u_prim_permuted[j, ii, 3]
      p_rr   = u_prim_permuted[j, ii, 4]

      # Compute required mean values
      rho_avg = 0.5 * (rho_ll + rho_rr)
      v1_avg  = 0.5 * ( v1_ll +  v1_rr)
      v2_avg  = 0.5 * ( v2_ll +  v2_rr)
      p_avg   = 0.5 * (  p_ll +   p_rr)
      kin_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)
      pv1_avg = 0.5 * (p_ll * v1_rr + p_rr * v1_ll)

      # Calculate fluxes depending on Cartesian orientation
      f1 = rho_avg * v1_avg
      f2 = f1 * v1_avg + p_avg
      f3 = f1 * v2_avg
      f4 = p_avg * v1_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv1_avg

      # Add scaled fluxes to RHS
      factor_i = alpha * derivative_split[i, ii]
      du_permuted[j, i, 1] += factor_i * f1
      du_permuted[j, i, 2] += factor_i * f2
      du_permuted[j, i, 3] += factor_i * f3
      du_permuted[j, i, 4] += factor_i * f4

      factor_ii = alpha * derivative_split[ii, i]
      du_permuted[j, ii, 1] += factor_ii * f1
      du_permuted[j, ii, 2] += factor_ii * f2
      du_permuted[j, ii, 3] += factor_ii * f3
      du_permuted[j, ii, 4] += factor_ii * f4
    end
  end

  @turbo for v in eachvariable(equations),
             j in eachnode(dg),
             i in eachnode(dg)
    du[i, j, v] = du_permuted[j, i, v]
  end

  # y direction
  # The memory layout is already optimal for SIMD vectorization in this loop.
  for j in eachnode(dg), jj in (j+1):nnodes(dg)
    @turbo for i in eachnode(dg)
      rho_ll = u_prim[i, j, 1]
      v1_ll  = u_prim[i, j, 2]
      v2_ll  = u_prim[i, j, 3]
      p_ll   = u_prim[i, j, 4]

      rho_rr = u_prim[i, jj, 1]
      v1_rr  = u_prim[i, jj, 2]
      v2_rr  = u_prim[i, jj, 3]
      p_rr   = u_prim[i, jj, 4]

      # Compute required mean values
      rho_avg = 0.5 * (rho_ll + rho_rr)
      v1_avg  = 0.5 * ( v1_ll +  v1_rr)
      v2_avg  = 0.5 * ( v2_ll +  v2_rr)
      p_avg   = 0.5 * (  p_ll +   p_rr)
      kin_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)
      pv2_avg = 0.5 * (p_ll * v2_rr + p_rr * v2_ll)

      # Calculate fluxes depending on Cartesian orientation
      f1 = rho_avg * v2_avg
      f2 = f1 * v1_avg
      f3 = f1 * v2_avg + p_avg
      f4 = p_avg*v2_avg * equations.inv_gamma_minus_one + f1 * kin_avg + pv2_avg

      # Add scaled fluxes to RHS
      factor_j = alpha * derivative_split[j, jj]
      du[i, j, 1] += factor_j * f1
      du[i, j, 2] += factor_j * f2
      du[i, j, 3] += factor_j * f3
      du[i, j, 4] += factor_j * f4

      factor_jj = alpha * derivative_split[jj, j]
      du[i, jj, 1] += factor_jj * f1
      du[i, jj, 2] += factor_jj * f2
      du[i, jj, 3] += factor_jj * f3
      du[i, jj, 4] += factor_jj * f4
    end
  end

  # Finally, we add the temporary RHS computed here to the global RHS in the
  # given `element`.
  @turbo for v in eachvariable(equations),
             j in eachnode(dg),
             i in eachnode(dg)
    _du[v, i, j, element] += du[i, j, v]
  end
end

@inline function split_form_kernel!(_du::PtrArray, u_cons::PtrArray,
                                    element, mesh::TreeMesh{2},
                                    nonconservative_terms::Val{false},
                                    equations::CompressibleEulerEquations2D,
                                    volume_flux::typeof(flux_ranocha_turbo),
                                    dg::DGSEM, cache, alpha)
  @unpack derivative_split = dg.basis

  # Create a temporary array that will be used to store the RHS with permuted
  # indices `[i, j, v]` to allow using SIMD instructions.
  # `StrideArray`s with purely static dimensions do not allocate on the heap.
  du = StrideArray{eltype(u_cons)}(undef,
    (ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))...,
     StaticInt(nvariables(equations))))

  # Convert conserved to primitive variables on the given `element`. In addition
  # to the usual primitive variables, we also compute logarithms of the density
  # and pressure to increase the performance of the required logarithmic mean
  # values.
  u_prim = StrideArray{eltype(u_cons)}(undef,
    (ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))...,
     StaticInt(nvariables(equations) + 2))) # We also compute "+ 2" logs

  @turbo for j in eachnode(dg), i in eachnode(dg)
    rho    = u_cons[1, i, j, element]
    rho_v1 = u_cons[2, i, j, element]
    rho_v2 = u_cons[3, i, j, element]
    rho_e  = u_cons[4, i, j, element]

    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))

    u_prim[i, j, 1] = rho
    u_prim[i, j, 2] = v1
    u_prim[i, j, 3] = v2
    u_prim[i, j, 4] = p
    u_prim[i, j, 5] = log(rho)
    u_prim[i, j, 6] = log(p)
  end

  # x direction
  # At first, we create new temporary arrays with permuted memory layout to
  # allow using SIMD instructions along the first dimension (which is contiguous
  # in memory).
  du_permuted = StrideArray{eltype(u_cons)}(undef,
    (StaticInt(nnodes(dg)), StaticInt(nnodes(dg)),
     StaticInt(nvariables(equations))))

  u_prim_permuted = StrideArray{eltype(u_cons)}(undef,
    (StaticInt(nnodes(dg)), StaticInt(nnodes(dg)),
     StaticInt(nvariables(equations) + 2)))

  @turbo for v in indices(u_prim, 3), # v in eachvariable(equations) misses +2 logs
             j in eachnode(dg),
             i in eachnode(dg)
    u_prim_permuted[j, i, v] = u_prim[i, j, v]
  end
  fill!(du_permuted, zero(eltype(du_permuted)))

  # Next, we basically inline the volume flux. To allow SIMD vectorization and
  # still use the symmetry of the volume flux and the derivative matrix, we
  # loop over the triangular part in an outer loop and use a plain inner loop.
  for i in eachnode(dg), ii in (i+1):nnodes(dg)
    @turbo for j in eachnode(dg)
      rho_ll     = u_prim_permuted[j, i, 1]
      v1_ll      = u_prim_permuted[j, i, 2]
      v2_ll      = u_prim_permuted[j, i, 3]
      p_ll       = u_prim_permuted[j, i, 4]
      log_rho_ll = u_prim_permuted[j, i, 5]
      log_p_ll   = u_prim_permuted[j, i, 6]

      rho_rr     = u_prim_permuted[j, ii, 1]
      v1_rr      = u_prim_permuted[j, ii, 2]
      v2_rr      = u_prim_permuted[j, ii, 3]
      p_rr       = u_prim_permuted[j, ii, 4]
      log_rho_rr = u_prim_permuted[j, ii, 5]
      log_p_rr   = u_prim_permuted[j, ii, 6]

      # Compute required mean values
      # We inline the logarithmic mean to allow LoopVectorization.jl to optimize
      # it efficiently. This is equivalent to
      #   rho_mean = ln_mean(rho_ll, rho_rr)
      x1 = rho_ll
      log_x1 = log_rho_ll
      y1 = rho_rr
      log_y1 = log_rho_rr
      x1_plus_y1 = x1 + y1
      y1_minus_x1 = y1 - x1
      z1 = y1_minus_x1^2 / x1_plus_y1^2
      special_path1 = x1_plus_y1 / (2 + z1*(2/3 + z1*(2/5 + 2/7*z1)))
      regular_path1 = y1_minus_x1 / (log_y1 - log_x1)
      rho_mean = ifelse(z1 < 1.0e-4, special_path1, regular_path1)

      # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
      # in exact arithmetic since
      #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
      #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
      # inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
      x2 = rho_ll * p_rr
      log_x2 = log_rho_ll + log_p_rr
      y2 = rho_rr * p_ll
      log_y2 = log_rho_rr + log_p_ll
      x2_plus_y2 = x2 + y2
      y2_minus_x2 = y2 - x2
      z2 = y2_minus_x2^2 / x2_plus_y2^2
      special_path2 = (2 + z2*(2/3 + z2*(2/5 + 2/7*z2))) / x2_plus_y2
      regular_path2 = (log_y2 - log_x2) / y2_minus_x2
      inv_rho_p_mean = p_ll * p_rr * ifelse(z2 < 1.0e-4, special_path2, regular_path2)

      v1_avg = 0.5 * (v1_ll + v1_rr)
      v2_avg = 0.5 * (v2_ll + v2_rr)
      p_avg  = 0.5 * ( p_ll +  p_rr)
      velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)

      # Calculate fluxes depending on Cartesian orientation
      f1 = rho_mean * v1_avg
      f2 = f1 * v1_avg + p_avg
      f3 = f1 * v2_avg
      f4 = f1 * ( velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one ) + 0.5 * (p_ll*v1_rr + p_rr*v1_ll)

      # Add scaled fluxes to RHS
      factor_i = alpha * derivative_split[i, ii]
      du_permuted[j, i, 1] += factor_i * f1
      du_permuted[j, i, 2] += factor_i * f2
      du_permuted[j, i, 3] += factor_i * f3
      du_permuted[j, i, 4] += factor_i * f4

      factor_ii = alpha * derivative_split[ii, i]
      du_permuted[j, ii, 1] += factor_ii * f1
      du_permuted[j, ii, 2] += factor_ii * f2
      du_permuted[j, ii, 3] += factor_ii * f3
      du_permuted[j, ii, 4] += factor_ii * f4
    end
  end

  @turbo for v in eachvariable(equations),
             j in eachnode(dg),
             i in eachnode(dg)
    du[i, j, v] = du_permuted[j, i, v]
  end

  # y direction
  # The memory layout is already optimal for SIMD vectorization in this loop.
  for j in eachnode(dg), jj in (j+1):nnodes(dg)
    @turbo for i in eachnode(dg)
      rho_ll     = u_prim[i, j, 1]
      v1_ll      = u_prim[i, j, 2]
      v2_ll      = u_prim[i, j, 3]
      p_ll       = u_prim[i, j, 4]
      log_rho_ll = u_prim[i, j, 5]
      log_p_ll   = u_prim[i, j, 6]

      rho_rr     = u_prim[i, jj, 1]
      v1_rr      = u_prim[i, jj, 2]
      v2_rr      = u_prim[i, jj, 3]
      p_rr       = u_prim[i, jj, 4]
      log_rho_rr = u_prim[i, jj, 5]
      log_p_rr   = u_prim[i, jj, 6]

      # Compute required mean values
      # We inline the logarithmic mean to allow LoopVectorization.jl to optimize
      # it efficiently. This is equivalent to
      #   rho_mean = ln_mean(rho_ll, rho_rr)
      x1 = rho_ll
      log_x1 = log_rho_ll
      y1 = rho_rr
      log_y1 = log_rho_rr
      x1_plus_y1 = x1 + y1
      y1_minus_x1 = y1 - x1
      z1 = y1_minus_x1^2 / x1_plus_y1^2
      special_path1 = x1_plus_y1 / (2 + z1*(2/3 + z1*(2/5 + 2/7*z1)))
      regular_path1 = y1_minus_x1 / (log_y1 - log_x1)
      rho_mean = ifelse(z1 < 1.0e-4, special_path1, regular_path1)

      # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
      # in exact arithmetic since
      #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
      #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
      # inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
      x2 = rho_ll * p_rr
      log_x2 = log_rho_ll + log_p_rr
      y2 = rho_rr * p_ll
      log_y2 = log_rho_rr + log_p_ll
      x2_plus_y2 = x2 + y2
      y2_minus_x2 = y2 - x2
      z2 = y2_minus_x2^2 / x2_plus_y2^2
      special_path2 = (2 + z2*(2/3 + z2*(2/5 + 2/7*z2))) / x2_plus_y2
      regular_path2 = (log_y2 - log_x2) / y2_minus_x2
      inv_rho_p_mean = p_ll * p_rr * ifelse(z2 < 1.0e-4, special_path2, regular_path2)

      v1_avg = 0.5 * (v1_ll + v1_rr)
      v2_avg = 0.5 * (v2_ll + v2_rr)
      p_avg  = 0.5 * ( p_ll +  p_rr)
      velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)

      # Calculate fluxes depending on Cartesian orientation
      f1 = rho_mean * v2_avg
      f2 = f1 * v1_avg
      f3 = f1 * v2_avg + p_avg
      f4 = f1 * ( velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one ) + 0.5 * (p_ll*v2_rr + p_rr*v2_ll)

      # Add scaled fluxes to RHS
      factor_j = alpha * derivative_split[j, jj]
      du[i, j, 1] += factor_j * f1
      du[i, j, 2] += factor_j * f2
      du[i, j, 3] += factor_j * f3
      du[i, j, 4] += factor_j * f4

      factor_jj = alpha * derivative_split[jj, j]
      du[i, jj, 1] += factor_jj * f1
      du[i, jj, 2] += factor_jj * f2
      du[i, jj, 3] += factor_jj * f3
      du[i, jj, 4] += factor_jj * f4
    end
  end

  # Finally, we add the temporary RHS computed here to the global RHS in the
  # given `element`.
  @turbo for v in eachvariable(equations),
             j in eachnode(dg),
             i in eachnode(dg)
    _du[v, i, j, element] += du[i, j, v]
  end
end


text:
module FTPServer

import Base: Process
import Base: close
using Conda
using Memento
using PyCall
using Random: randstring
const LOGGER = getlogger(@__MODULE__)

const pyopenssl_crypto = PyNULL()
const pyopenssl_SSL = PyNULL()
const pyftpdlib_authorizers = PyNULL()
const pyftpdlib_handlers = PyNULL()
const pyftpdlib_servers = PyNULL()

# Defaults from pyftpdlib example
const HOST = "localhost"
const PERM = "elradfmwM"

const SCRIPT = abspath(dirname(@__FILE__), "server.py")
const ROOT = abspath(joinpath(dirname(dirname(@__FILE__)), "deps", "usr", "ftp"))
const HOMEDIR = joinpath(ROOT, "data")
const CERT = joinpath(ROOT, "test.crt")
const KEY = joinpath(ROOT, "test.key")

function __init__()
    Memento.register(LOGGER)

    copy!(pyopenssl_crypto, pyimport_conda("OpenSSL.crypto", "OpenSSL"))
    copy!(pyopenssl_SSL, pyimport_conda("OpenSSL.SSL", "OpenSSL"))

    # Note: For `pyftpdlib` we'll specify an exact version to make behaviour of FTPServer.jl
    # consistent when rolling back to an earlier version.
    # For details see: https://github.com/invenia/FTPClient.jl/issues/91#issuecomment-632698841
    copy!(pyftpdlib_servers, pyimport_conda("pyftpdlib.servers", "pyftpdlib==1.5.6", "invenia"))

    mkpath(ROOT)
end

mutable struct Server
    homedir::AbstractString
    port::Int
    username::AbstractString
    password::AbstractString
    permissions::AbstractString
    security::Symbol
    process::Process
    io::IO
end

"""
    Server(
        homedir::AbstractString=$HOMEDIR;
        username::AbstractString="",
        password::AbstractString="",
        permissions::AbstractString=$PERM,
        security::Symbol=:none,
        force_gen_certs::Bool=true,
        debug_command::Bool=false,
    )

A Server stores settings to create a pyftpdlib server.

# Arguments
- `homedir::AbstractString=$HOMEDIR`: Directory where you want to store your data for the
  test server.

# Keywords
- `username::AbstractString=""`: Default login username. Defaults to 'userXXXX' where 'XXXX'
  is a number between 1 and 9999.
- `password::AbstractStringi=""`: Default login password. Defalts to a random string of 40
  characters.
- `permission::AbstractString=$PERM`: Default user read/write permissions.
- `security::Symbol=:none`: Security method to use for connecting (options: `:none`,
  `:implicit`, `:explicit`). Passing in `:none` will use FTP and passing in `:implicit` or
  `:explicit` will use the appropriate FTPS connection.
- `force_gen_certs::Bool=true`: Force regeneration of certificate and key file.
- `debug_command::Bool=false`: Print out the python command being used, for debugging
  purposes.
"""
function Server(
    homedir::AbstractString=HOMEDIR;
    username::AbstractString="",
    password::AbstractString="",
    permissions::AbstractString=PERM,
    security::Symbol=:none,
    force_gen_certs::Bool=true,
    debug_command::Bool=false,
)
    if isempty(username)
        username = string("user", rand(1:9999))
    end
    if isempty(password)
        password = randstring(40)
    end

    python = PyCall.current_python()
    cmd = `$python $SCRIPT $username $password $homedir --permissions $permissions`
    if security != :none
        cmd = `$cmd --tls $security --cert-file $CERT --key-file $KEY`

        if force_gen_certs
            cmd = `$cmd --force-gen-certs`
        end
    end

    # If we're having issues with the above command, it can be useful to print it out so we
    # can run the command against the python script itself to see if it gives us any extra
    # insight.
    if debug_command
        info(LOGGER, "Running Command: $cmd")
    end

    # Note: open(::AbstractCmd, ...) won't work here as it doesn't allow us to capture
    # STDERR.
    io = Pipe()
    process = run(pipeline(cmd, stdout=io, stderr=io), wait=false)

    # Grab the Port value from the python script output
    line = readline(io)
    while !occursin("starting FTP", line)
        line = readline(io)
    end
    m = match(r"starting FTP.* server on .*:(?<port>\d+)", line)

    # If we found the port, store the server data in an object, else show an error
    if m !== nothing
        port = parse(Int, m[:port])
        Server(homedir, port, username, password, permissions, security, process, io)
    else
        kill(process)
        error(line, String(readavailable(io)))  # Display traceback
    end
end

"""
    serve(f, args...; kwargs...)

Passes `args` and `kwargs` to the `Server` constructor and runs the function `f` by passing
in the `server` instance. Upon completion the `server` will automatically be shutdown.
"""
function serve(f::Function, args...; kwargs...)
    server = Server(args...; kwargs...)

    try
        f(server)
    finally
        close(server)
    end
end

hostname(server::Server) = HOST
port(server::Server) = server.port
username(server::Server) = server.username
password(server::Server) = server.password
close(server::Server) = kill(server.process)

localpath(server::Server, path::AbstractString) = joinpath(
    server.homedir, split(path, '/')...
)

function tempfile(path::AbstractString)
    content = randstring(rand(1:100))
    open(path, "w") do fp
        write(fp, content)
    end
    return content
end

function setup_home(dir::AbstractString)
    mkdir(dir)
    tempfile(joinpath(dir, "test_download.txt"))
    tempfile(joinpath(dir, "test_download2.txt"))
    mkdir(joinpath(dir, "test_directory"))
end

"""
    init()

Creates a test $HOMEDIR with a few sample files if one hasn't already been setup.

```
$HOMEDIR/test_download.txt
$HOMEDIR/test_download2.txt
$HOMEDIR/test_directory/
````
"""
init() = isdir(FTPServer.HOMEDIR) || setup_home(FTPServer.HOMEDIR)

"""
    cleanup()

Cleans up the default FTPServer.ROOT directory:

- $HOMEDIR
- $CERT
- $KEY
"""
function cleanup()
    rm(FTPServer.HOMEDIR, recursive=true)
    isfile(FTPServer.CERT) && rm(FTPServer.CERT)
    isfile(FTPServer.KEY) && rm(FTPServer.KEY)
end

"""
    uri(server::Server)

Create an FTP URI from an FTP server object.

# Arguments
- `server::Server`: FTPServer object
"""
function uri(server::Server)
    ftp_type = if server.security === :implicit
        "ftps"
    server.security === :explicit
        "ftpes"
    else
        "ftp"
    end

    string(
        "$ftp_type://$(username(server)):$(password(server))",
        "@$(hostname(server)):$(port(server))",
    )
end

end # module


text:
struct GradientDescent{T1} <: QuasiNewton{T1}
    approx::T1
end
update!(scheme::GradientDescent{<:Inverse}, A::UniformScaling, s, y) = A
update(scheme::GradientDescent{<:Inverse}, A::UniformScaling, s, y) = A
update!(scheme::GradientDescent{<:Direct}, A, s, y) = A
update(scheme::GradientDescent{<:Direct}, A, s, y) = A


text:
export State, KeplerianState, array, epoch, scale, frame, body, epoch_type, pos_type, vel_type

#################
# AbstractState #
#################

export AbstractState, epoch, frame

abstract type AbstractState{Scale, Frame, Body} end

function epoch end

scale(s::AbstractState{Scale}) where {Scale} = Scale
frame(s::AbstractState{_S, Frame}) where {_S, Frame} = Frame
body(s::AbstractState{_S, _F, Body}) where {_S, _F, Body} = Body
array(s::AbstractState) = vcat(state(s)...)

function period(s::AbstractState)
    μ = grav_param(body(s))
    a, _ = keplerian(s)
    return period(a, μ)
end

struct State{Scale, Frame, Body, T, TP, TV} <: AbstractState{Scale, Frame, Body}
    epoch::Epoch{Scale, T}
    pos::SVector{3, TP}
    vel::SVector{3, TV}
    function State(epoch::Epoch{Scale, T}, pos, vel;
                   frame::Frame=icrf, body::Body=earth) where {Scale, Frame, Body, T}
        TP = eltype(pos)
        TV = eltype(vel)
        new{Scale::TimeScale, frame::AbstractFrame, body::CelestialBody, T, TP, TV}(epoch, pos, vel)
    end
end

function (==)(s1::State, s2::State)
    s1.epoch == s2.epoch && s1.pos == s2.pos && s1.vel == s2.vel
end

function Base.isapprox(s1::State, s2::State)
    s1.epoch ≈ s2.epoch && s1.pos ≈ s2.pos && s1.vel ≈ s2.vel
end

epoch_type(s::State{_S, _F, _B, T}) where {_S, _F, _B, T} = T
pos_type(s::State{_S, _F, _B, _T, TP}) where {_S, _F, _B, _T, TP} = TP
vel_type(s::State{_S, _F, _B, _T, _TP, TV}) where {_S, _F, _B, _T, _TP, TV} = TV

epoch(s::State) = s.epoch
position(s::State) = s.pos
velocity(s::State) = s.vel
state(s::State) = (s.pos, s.vel)
keplerian(s::State) = keplerian(position(s), velocity(s), grav_param(body(s)))

struct KeplerianState{Scale, Frame, Body, T} <: AbstractState{Scale, Frame, Body}
    epoch::Epoch{Scale, T}
    a::T
    e::T
    i::T
    Ω::T
    ω::T
    ν::T
    function KeplerianState(epoch::Epoch{Scale, T}, a, e, i, Ω, ω, ν;
                            frame::Frame=icrf, body::Body=earth) where {Scale, Frame, Body, T}
        new{Scale::TimeScale, frame::AbstractFrame, body::CelestialBody, T}(epoch, a, e, i, Ω, ω, ν)
    end
end

function (==)(s1::KeplerianState, s2::KeplerianState)
    return s1.epoch == s2.epoch &&
        s1.a == s2.a &&
        s1.e == s2.e &&
        s1.i == s2.i &&
        s1.Ω == s2.Ω &&
        s1.ω == s2.ω &&
        s1.ν == s2.ν
end

function Base.isapprox(s1::KeplerianState, s2::KeplerianState)
    return s1.epoch ≈ s2.epoch &&
        s1.a ≈ s2.a &&
        s1.e ≈ s2.e &&
        s1.i ≈ s2.i &&
        s1.Ω ≈ s2.Ω &&
        s1.ω ≈ s2.ω &&
        s1.ν ≈ s2.ν
end

Base.eltype(::KeplerianState{_S, _F, _B, T}) where {_S, _F, _B, T} = T

epoch(s::KeplerianState) = s.epoch
position(s::KeplerianState) = cartesian(s.a, s.e, s.i, s.Ω, s.ω, s.ν, grav_param(body(s)))[1]
velocity(s::KeplerianState) = cartesian(s.a, s.e, s.i, s.Ω, s.ω, s.ν, grav_param(body(s)))[2]
state(s::KeplerianState) = cartesian(s.a, s.e, s.i, s.Ω, s.ω, s.ν, grav_param(body(s)))
keplerian(s::KeplerianState) = (s.a, s.e, s.i, s.Ω, s.ω, s.ν)

KeplerianState(s::AbstractState) = KeplerianState(epoch(s), keplerian(s)...; frame=frame(s), body=body(s))
State(s::KeplerianState) = State(epoch(s), state(s)...; frame=frame(s), body=body(s))

# struct ThreeBodyState{
#         F<:Frame,
#         T,
#         C1<:CelestialBody,
#         C2<:CelestialBody,
#     } <: AbstractState
#     ep::Epoch{T}
#     r::SVector{3,Float64}
#     v::SVector{3,Float64}
#
#     function ThreeBodyState(
#         ep::Epoch{T}, r, v, frame::Type{F}=GCRF,
#         primary::Type{C1}=Sun, secondary::Type{C2}=Earth
#     ) where {T,F<:Frame,C1<:CelestialBody,C2<:CelestialBody}
#         new{F, T, C1, C2}(ep, r, v)
#     end
# end
#
# frame(::ThreeBodyState{F}) where F<:Frame = F
# timescale(::ThreeBodyState{<:Frame, T}) where T = T
# primary(::ThreeBodyState{<:Frame,,C}) where C<:CelestialBody = C
# secondary(::ThreeBodyState{<:Frame,,<:CelestialBody,C}) where {
#     C<:CelestialBody} = C


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using Test

using QuantumCircuits
using QuantumCircuits.QCircuits.Gates
using QuantumCircuits.QCircuits.Math
using QuantumCircuits.QCircuits.Circuit: toQiskit, getCode

################################################################################
#  Simplify - one gate                                                         #
################################################################################
const u2(qubit, ϕ, λ) = (U3, (qubit, π/2, ϕ, λ))
const u1(qubit, λ) = (U3, (qubit, 0, 0, λ))

function test_simplify(u, args, gate)
    qc = QCircuit(1)
    add!(qc, u, args...)
    sqc = simplify(qc)
    code = getCode(sqc)
    @test length(code) == 1
    @test typeof(code[1]) == gate
end
test_simplify(U3, (0, π, 0, π), X)
test_simplify(U3, (0, π, π/2, π/2), Y)
test_simplify(u1(0, π)..., Z)
test_simplify(u2(0, 0, π)..., H)
test_simplify(u1(0, π/2)..., S)
test_simplify(u1(0, -π/2)..., Sd)
test_simplify(u1(0, π/4)..., T)
test_simplify(u1(0, -π/4)..., Td)

test_simplify(U3, (0, -π, 3.1415, 0.004643), X)
# qc = QCircuit(1)
# qc.u3(0, -π, 3.1415, 0.004643)
# println(simplify(qc))
# using QGen.Quantum.Gates: Xmatrix
# unitary_error(tomatrix(qc), Xmatrix)

################################################################################
#  Simplify - identify                                                         #
################################################################################
qc = QCircuit(1)
qc.u3(0, 0, -3.678862306482877, -2.6043230006967093)
sqc = simplify(qc)
@test length(getCode(sqc)) == 0

################################################################################
#  Simplify - 2 gates                                                          #
################################################################################
function test_simplify_gates(u3, args, gates, orGate=nothing)
    qc = QCircuit(1)
    add!(qc, u3, args...)
    sqc = simplify(qc)
    code = getCode(sqc)
    @test length(code) == length(gates)
    for (i, g) in enumerate(gates)
        if isnothing(orGate)
            @test typeof(code[i]) == g
        else
            @test typeof(code[i]) == g || typeof(code[i]) == orGate[i]
        end
    end

    @test unitary_error(tomatrix(qc), tomatrix(sqc)) < 1e-8
end

test_simplify_gates(U3, (0, -π/2, π, π), [H, X], [Z, H])
test_simplify_gates(U3, (0, π, 0, π/2), [Y, Sd], [Sd, X])
test_simplify_gates(U3, (0, π, 0, 5π/4), [T, X])
test_simplify_gates(U3, (0, π, 2π, 0), [Z, X])
test_simplify_gates(U3, (0, 0, 2.3114387473071725, 2.400950233077517), [Sd])

# qc = QCircuit(1)
# add!(qc, U3(0, 0, 2.3114387473071725, 2.400950233077517))
# tomatrix(qc)
#
#
# sqc = simplify(qc)
#
# qc1 = QCircuit(1)
# add!(qc1, Z(0))
# add!(qc1, S(0))
# tomatrix(qc1)
# unitary_error(tomatrix(qc), tomatrix(qc1))


text:
module TestMathProgBase

import Mosek.MosekMathProgSolverInterface
import Mosek
using MathProgBase
using Base.Test

#const MPB = MathProgBase
include(joinpath(Pkg.dir("MathProgBase","test","linprog.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","mixintprog.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","quadprog.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","nlp.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","conicinterface.jl")))
include(joinpath(Pkg.dir("MathProgBase","test","linproginterface.jl")))

function test_mathprogbase(solver)
    @testset "linprog" begin
        linprogtest(solver)
    end

    @testset "mip" begin
        mixintprogtest(solver)
    end

    @testset "qp" begin
        #quadprogtest(solver)
        #qpdualtest(solver)
        #socptest(solver)
    end

    @testset "nlp" begin
        #nlptest(solver)
        #nlptest_nohessian(solver)
        #convexnlptest(solver)
        #rosenbrocktest(solver)
    end

    @testset "conic" begin
        coniclineartest(solver, duals=true)
        conicSOCtest(solver, duals=true)
        conicSOCRotatedtest(solver, duals=true)
        conicSOCINTtest(solver)
        conicEXPtest(solver, duals=true)
        conicSDPtest(solver, duals=true)
    end

    @testset "linproginterface" begin
        linprogsolvertest(solver)
        linprogsolvertestextra(solver)
    end
end

test_mathprogbase(Mosek.MosekSolver(QUIET = true))

end


text:
"""Only used when building the model"""
function _retrieve_forecasts(sys::PSY.System, ::Type{C}) where {C <: PSY.Component}

    first_step = PSY.get_forecasts_initial_time(sys)
    forecasts = PSY.get_component_forecasts(C, sys, first_step)
    isempty(forecasts) && error("System has no forecasts for device $(C)")

    return collect(forecasts)
    #TODO: Remove this collect
end


text:
using BenchmarkTools
const bench = SUITE = BenchmarkGroup()

using Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.WiringDiagrams
using Catlab.WiringDiagrams.DirectedWiringDiagrams: TheoryWiringDiagram
using Catlab.Graphs.BasicGraphs: TheoryGraph
using Random: rand
using CSetAutomorphisms

# Color refine benchmarks
#########################

@acset_type WD(TheoryWiringDiagram)

# WDs
#----

# This could be generalized to a function making a 'random' C-set.
sizes = (Box=3, Wire=3, InPort=3, OutPort=3, OuterInPort=3, OuterOutPort=3,
         OutWire=3, InWire=3, PassWire=3)
wd = WD()
for (k,v) in pairs(sizes)
  add_parts!(wd, k, v)
end
for g in TheoryWiringDiagram.generators[:Hom]
  n, d_, cd_ = g.args
  d, cd = d_.args[1], cd_.args[1]
  set_subpart!(wd, n, rand(1:sizes[cd], sizes[d]))
end
[copy_parts!(wd, wd) for _ in 1:3] # add 8-fold symmetry

bench["Color saturate a random wiring diagram:"] = @benchmarkable color_saturate(wd)

# Graphs
#-------
g1 = star_graph(Graph, 4)
g2 = path_graph(Graph, 4)
g = Graph(1)
[copy_parts!(g, h) for h in [g1, g2, g1, g2]]

bench["Color saturate a graph:"]= @benchmarkable color_saturate(g)

G = Graph(1)
[copy_parts!(G, x) for x in [g1,g2]]
canonical_hash(G); canonical_hash(G; pres=TheoryGraph)
bench["Canonical form, nauty:"]= @benchmarkable canonical_hash(G; pres=TheoryGraph)
bench["Canonical form, csets:"]= @benchmarkable canonical_hash(G)


text:
include("gemv.jl")
include("gbmv.jl")
include("hemv.jl")
include("hbmv.jl")
include("symv.jl")
include("sbmv.jl")
include("trmv.jl")
include("trsv.jl")
include("ger.jl")
include("her.jl")
include("syr.jl")


text:
using Test
using Bigsimr
using Polynomials
using Distributions

@testset "Pearson Correlation Matching" begin
    dA = Beta(2, 3)
    dB = Binomial(2, 0.2)
    dC = Binomial(20, 0.2)

    @testset "Continuous-Continuous" begin
        @test -0.914 ≈ pearson_match(-0.9, dA, dA) atol=0.05
        @test -0.611 ≈ pearson_match(-0.6, dA, dA) atol=0.05
        @test -0.306 ≈ pearson_match(-0.3, dA, dA) atol=0.05
        @test  0.304 ≈ pearson_match( 0.3, dA, dA) atol=0.05
        @test  0.606 ≈ pearson_match( 0.6, dA, dA) atol=0.05
        @test  0.904 ≈ pearson_match( 0.9, dA, dA) atol=0.05
    end

    @testset "Discrete-Discrete" begin
        @test -0.937 ≈ pearson_match(-0.5, dB, dB) atol=0.05
        @test -0.501 ≈ pearson_match(-0.3, dB, dB) atol=0.05
        @test -0.322 ≈ pearson_match(-0.2, dB, dB) atol=0.05
        @test  0.418 ≈ pearson_match( 0.3, dB, dB) atol=0.05
        @test  0.769 ≈ pearson_match( 0.6, dB, dB) atol=0.05
        @test  0.944 ≈ pearson_match( 0.8, dB, dB) atol=0.05

        @test -0.939 ≈ pearson_match(-0.9, dC, dC) atol=0.05
        @test -0.624 ≈ pearson_match(-0.6, dC, dC) atol=0.05
        @test -0.311 ≈ pearson_match(-0.3, dC, dC) atol=0.05
        @test  0.310 ≈ pearson_match( 0.3, dC, dC) atol=0.05
        @test  0.618 ≈ pearson_match( 0.6, dC, dC) atol=0.05
        @test  0.925 ≈ pearson_match( 0.9, dC, dC) atol=0.05
    end

    @testset "Mixed" begin
        @test -0.890 ≈ pearson_match(-0.7, dB, dA) atol=0.05
        @test -0.632 ≈ pearson_match(-0.5, dB, dA) atol=0.05
        @test -0.377 ≈ pearson_match(-0.3, dB, dA) atol=0.05
        @test  0.366 ≈ pearson_match( 0.3, dB, dA) atol=0.05
        @test  0.603 ≈ pearson_match( 0.5, dB, dA) atol=0.05
        @test  0.945 ≈ pearson_match( 0.8, dB, dA) atol=0.05

        @test -0.928 ≈ pearson_match(-0.9, dC, dA) atol=0.05
        @test -0.618 ≈ pearson_match(-0.6, dC, dA) atol=0.05
        @test -0.309 ≈ pearson_match(-0.3, dC, dA) atol=0.05
        @test  0.308 ≈ pearson_match( 0.3, dC, dA) atol=0.05
        @test  0.613 ≈ pearson_match( 0.6, dC, dA) atol=0.05
        @test  0.916 ≈ pearson_match( 0.9, dC, dA) atol=0.05
    end

end

@testset "Pearson Correlation Utilities" begin

    @testset "Get Hermite Coefficients" begin
        dA = Binomial(20, 0.2)
        dB = NegativeBinomial(20, 0.002)
        dC = LogitNormal(3, 1)
        dD = Beta(5, 3)

        @test_nowarn Bigsimr.get_coefs(dA, 7)
        @test_nowarn Bigsimr.get_coefs(dB, 7)
        @test_nowarn Bigsimr.get_coefs(dC, 7)
        @test_nowarn Bigsimr.get_coefs(dD, 7)

        @test_nowarn Bigsimr.get_coefs(dA, 7.0)
        @test_nowarn Bigsimr.get_coefs(dB, 7.0)
        @test_nowarn Bigsimr.get_coefs(dC, 7.0)
        @test_nowarn Bigsimr.get_coefs(dD, 7.0)

        @test_throws InexactError Bigsimr.get_coefs(dA, 7.5)
        @test_throws InexactError Bigsimr.get_coefs(dB, 7.5)
        @test_throws InexactError Bigsimr.get_coefs(dC, 7.5)
        @test_throws InexactError Bigsimr.get_coefs(dD, 7.5)
    end

    @testset "Core Hermite Function" begin
        # Must work for any real input
        test_types = (Float64, Float32, Float16, BigFloat, Int128, Int64, Int32, Int16, BigInt, Rational)
        for T in test_types
            @test_nowarn Bigsimr._h(one(T), 5)
        end

        # For the following types, the input type should be the same as the output
        test_types = (Float64, Float32, Float16, BigFloat, Int64, BigInt, Rational)
        for T in test_types
            @test Bigsimr._h(one(T), 5) isa T
        end

        @test_nowarn Bigsimr._h(3.14, 5.0)
        @test_throws InexactError Bigsimr._h(3.14, 5.5)

        # Must work for arrays/matrices/vectors
        A = rand(3)
        B = rand(3, 3)
        C = rand(3, 3, 3)
        @test_nowarn Bigsimr._h(A, 3)
        @test_nowarn Bigsimr._h(B, 3)
        @test_nowarn Bigsimr._h(C, 3)
    end

    @testset "Hermite-Normal PDF" begin
        @test iszero(Bigsimr.Hp(Inf, 10))
        @test iszero(Bigsimr.Hp(-Inf, 10))
        @test 1.45182435 ≈ Bigsimr.Hp(1.0, 5)
    end

    @testset "Solve Polynomial on [-1, 1]" begin
        r1 = -1.0
        r2 = 1.0
        r3 = eps()
        r4 = 2 * rand() - 1

        P1 = coeffs(3 * fromroots([r1, 7, 7, 8]))
        P2 = coeffs(-5 * fromroots([r2, -1.14, -1.14, -1.14, -1.14, 1119]))
        P3 = coeffs(1.2 * fromroots([r3, nextfloat(1.0), prevfloat(-1.0)]))
        P4 = coeffs(fromroots([-5, 5, r4]))
        P5 = coeffs(fromroots([nextfloat(1.0), prevfloat(-1.0)]))
        P6 = coeffs(fromroots([-0.5, 0.5]))

        # One root at -1.0
        @test Bigsimr.solve_poly_pm_one(P1) ≈ r1 atol=0.001
        # One root at 1.0
        @test Bigsimr.solve_poly_pm_one(P2) ≈ r2 atol=0.001
        # Roots that are just outside [-1, 1]
        @test Bigsimr.solve_poly_pm_one(P3) ≈ r3 atol=0.001
        @test Bigsimr.solve_poly_pm_one(P4) ≈ r4 atol=0.001
        # Case of no roots
        @test isnan(Bigsimr.solve_poly_pm_one(P5))
        # Case of multiple roots
        @test length(Bigsimr.solve_poly_pm_one(P6)) == 2
        @test Bigsimr.nearest_root(-0.6, Bigsimr.solve_poly_pm_one(P6)) ≈ -0.5 atol=0.001
    end

end


text:
function raindrops(number)
    # has 3 as a factor, add 'Pling' to the result.
    # has 5 as a factor, add 'Plang' to the result.
    # has 7 as a factor, add 'Plong' to the result.
    # does not have any of 3, 5, or 7 as a factor, the result should be the digits of the number.
    result = ""
    for (factor, string) in [(3, "Pling"), (5, "Plang"), (7, "Plong")]
        if number % factor == 0
            result *= string
        end
    end
    if result == ""
        return "$number"
    else
        return result
    end
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
#=
MNIST sample
Taken from model zoo of Julia and modified.
https://github.com/FluxML/model-zoo/blob/master/vision/mnist/mlp.jl
=#

using Printf

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, @epochs
using Base.Iterators: partition
using BSON: @load, @save
using CuArrays
using Random

function prepare_dataset(;train=true)
    train_or_test = ifelse(train,:train,:test)
    imgs = MNIST.images(train_or_test)
    X = hcat(float.(vec.(imgs))...)
    labels = MNIST.labels(train_or_test)
    Y = onehotbatch(labels,0:9)
    return X, Y
end

function define_model(;hidden)
    mlp = Chain(Dense(28^2,hidden,relu),
                Dense(hidden,hidden,relu),
                Dense(hidden,10),
                softmax)
    return mlp
end

function split_dataset_random(X, Y)
    divide_ratio=0.9
    shuffled_indices = shuffle(1:size(Y)[2])
    divide_idx = round(Int,0.9*length(shuffled_indices))
    train_indices = shuffled_indices[1:divide_idx]
    val_indices = shuffled_indices[divide_idx:end]
    train_X = X[:,train_indices]
    train_Y = Y[:,train_indices]
    val_X = X[:,val_indices]
    val_Y = Y[:,val_indices]
    return train_X, train_Y, val_X, val_Y
end

function train()
    println("Start to train")
    epochs = 10
    X, Y = prepare_dataset(train=true)
    train_X, train_Y, val_X,val_Y = split_dataset_random(X, Y)
    model = define_model(hidden=100) |> gpu
    loss(x,y)= crossentropy(model(x),y)
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    batchsize = 64
    train_dataset = gpu.([(train_X[:,batch] ,train_Y[:,batch]) for batch in partition(1:size(train_Y)[2],batchsize)])
    val_dataset = gpu.([(val_X[:,batch] ,val_Y[:,batch]) for batch in partition(1:size(val_Y)[2],batchsize)])

    callback_count = 0
    eval_callback = function callback()
        callback_count += 1
        if callback_count == length(train_dataset)
            println("action for each epoch")
            total_loss = 0
            total_acc = 0
            for (vx, vy) in val_dataset
                total_loss += loss(vx, vy)
                total_acc += accuracy(vx, vy)
            end
            total_loss /= length(val_dataset)
            total_acc /= length(val_dataset)
            @show total_loss, total_acc
            callback_count = 0
            pretrained = model |> cpu
            @save "pretrained.bson" pretrained
            callback_count = 0
        end
        if callback_count % 50 == 0
            progress = callback_count / length(train_dataset)
           @printf("%.3f\n", progress)
        end
    end
    optimizer = ADAM(params(model))

    @epochs epochs Flux.train!(loss, train_dataset, optimizer, cb = eval_callback)

    pretrained = model |> cpu
    weights = Tracker.data.(params(pretrained))
    @save "pretrained.bson" pretrained
    @save "weights.bson" weights
    println("Finished to train")
end

function predict()
    println("Start to evaluate testset")
    println("loading pretrained model")
    @load "pretrained.bson" pretrained
    model = pretrained |> gpu
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    println("prepare dataset")
    X, Y = prepare_dataset(train=false)
    X = X |> gpu
    Y = Y |> gpu
    @show accuracy(X, Y)
    println("Done")
end

function predict2()
    println("Start to evaluate testset")
    println("loading pretrained model")
    @load "weights.bson" weights
    model = define_model(hidden=100)
    Flux.loadparams!(model, weights)
    model = model |> gpu
    accuracy(x, y) = mean(onecold(model(x)) .== onecold(y))
    println("prepare dataset")
    X, Y = prepare_dataset(train=false)
    X = X |> gpu
    Y = Y |> gpu
    @show accuracy(X, Y)
    println("Done")
end

function main()
    train()
    predict()
    predict2()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
#==============================================================================
    colorscheme: cathode_basic
        REPL colorschme definition
        load this during julia startup (not extending OhMyREPL)
        iterm color: Cathode*

    usage:
                include("/Users/merlinr/repo/julia-pkgs/CathodeJulia/startup-jl-colorscheme-cathode_basic.jl")
    tested:     NO
    author:     merlinr@
 =============================================================================#

using Crayons

using OhMyREPL: Passes.SyntaxHighlighter

function create_cathode_basic()
    cs = SyntaxHighlighter.ColorScheme()
    SyntaxHighlighter.symbol!(cs, Crayon(foreground = :light_red, bold = true))
    SyntaxHighlighter.macro!(cs, Crayon(foreground = :light_cyan))
    SyntaxHighlighter.keyword!(cs, Crayon(foreground = :cyan, bold = true))
    SyntaxHighlighter.function_def!(cs, Crayon(foreground = :magenta, bold = true))
    SyntaxHighlighter.call!(cs, Crayon(foreground = :light_magenta))
    SyntaxHighlighter.op!(cs, Crayon(foreground = :light_blue))
    SyntaxHighlighter.number!(cs, Crayon(foreground = :light_yellow))
    SyntaxHighlighter.text!(cs, Crayon(foreground = :green))
    SyntaxHighlighter.comment!(cs, Crayon(foreground = :black))
    SyntaxHighlighter.string!(cs, Crayon(foreground = :red))
    SyntaxHighlighter.argdef!(cs, Crayon(foreground = :blue))
    SyntaxHighlighter.error!(cs, Crayon(foreground = :light_gray))
    return cs
end

SyntaxHighlighter.add!("CathodeBasic", create_cathode_basic())

colorscheme!("CathodeBasic")

# remove the rainbow
# OhMyREPL.enable_pass!("RainbowBrackets", false)

# display the REPL state
OhMyREPL.showpasses()

max_stars_repo_name: mahiki/CathodeJulia

max_issues_repo_name: mahiki/CathodeJulia

max_forks_repo_name: mahiki/CathodeJulia


text:
function fib_producer(c::Channel)
    a, b = (0, 1)
    for i = 1:10
        put!(c, b)
        a, b = (b, a + b)
    end
end

chnl = Channel(fib_producer) #> Channel{Any}(sz_max:0,sz_curr:1)

take!(chnl) #> 1
take!(chnl) #> 1
take!(chnl) #> 2
take!(chnl) #> 3
take!(chnl) #> 5
take!(chnl) #> 8
take!(chnl) #> 13
take!(chnl) #> 21
take!(chnl) #> 34
take!(chnl) #> 55
# take!(chnl) #> ERROR: InvalidStateException("Channel is closed.", :closed)

for n in chnl
    println(n)
end
#> 1 1 1 2 3 5 8 13 21 34 55

chnl = @task fib_producer(c::Channel) # Task (runnable) @0x0000000005696d80

fac(i::Integer) = (i > 1) ? i*fac(i - 1) : 1
c = Channel(0)
task = @async foreach(i->put!(c,fac(i)), 1:5)
bind(c,task)
for i in c
   @show i
end
# i = 1
# i = 2
# i = 6
# i = 24
# i = 120


text:
using Dates
using PcapTools
using Test
using UnixTimes

@testset "PcapTools.jl" begin
    include("reader_tests.jl")
    include("writer_tests.jl")
end


text:
"""

Package `FITSexplore`

"""
module FITSexplore

export fitsexplore

using FITSIO, EasyFITS, ArgParse

const suffixes = [".fits", ".fits.gz","fits.Z",".oifits",".oifits.gz",".oifits.Z"]

"""
	endswith(chain::Union{String,Vector{String}}, pattern::Vector{String})

Overloading of Base.endswith for vectors of string. Return `true` if `chain`
ends with one of the patterns given in `pattern`
"""
function Base.endswith(chain::Union{String,Vector{String}}, pattern::Vector{String})
	for str in pattern
		if endswith(chain,str)
			return true
		end
	end
	return false
end

"""
	endswith(chain::Vector{String}, pattern::String)

Overloading of Base.endswith for vectors of string. Return `true` if `chain`
ends with `pattern`
"""
function Base.endswith(chains::Vector{String},pattern::AbstractString)
	for chain in chains
		if endswith(chain,pattern)
			return true
		end
	end
	return false
end

"""
	newlist = filtercat(filelist,keyword,value)

Build a `newlist` dictionnary of all files where `fitsheader[keyword] == value`.
"""
function filtercat(filelist::Dict{String, FitsHeader},
					keyword::String,
					values::Union{Vector{String}, Vector{Bool}, Vector{Integer}, Vector{AbstractFloat}})
	newlist = Dict{String, FitsHeader}()
	for value in values
		merge!(newlist, filtercat(filelist,keyword,value))
	end
	return newlist
end

function filtercat(filelist::Dict{String, FitsHeader},
					keyword::String,
					value::Union{String, Bool, Integer, AbstractFloat, Nothing})
	try tmp = filter(p->p.second[keyword] == value,filelist)
		return  tmp
	catch
		return Dict{String, FitsHeader}()
	end
end

function fitsexplore(dir::String)
	filedict = Dict{String, FitsHeader}()

	for filename in readdir(dir,join=true)
		if isfile(filename)
			if endswith(filename,suffixes)
				get!(filedict, filename) do
					read(FitsHeader, filename)
				end
			end
		end
	end
	return filedict
end

function parse_keywords(args::Vector{String}, keywords::Vector{Vector{String}} )
	parse_keywords(args,[k[1] for k in keywords])
end

function parse_keywords(args::Vector{String}, keywords::Vector{String} )
	for filename in args
		if isfile(filename)
			if endswith(filename,suffixes)
				header  = read(FitsHeader, filename)
				str = ""
				iskeyword = true
				for key in keywords
					if haskey(header,key)
						str = str * "\t" * string(header[key])
					else
						iskeyword = false
					end
				end
				if iskeyword
					println(filename ,"\t", str)
				end
			end
		end
	end
end

function parse_filter(args::Vector{String}, filter::Vector{String} )
	for filename in args
		if isfile(filename)
			if endswith(filename,suffixes)
				header  = read(FitsHeader, filename)
				if haskey(header,filter[1])
					if comparekeys(header[filter[1]],filter[2])
						println(filename)
					end
				end
			end
		end
	end
end

function comparekeys(key1::AbstractString,key2::AbstractString)
	return key1==key2
end

function comparekeys(key1::Bool,key2::AbstractString)
	p = false
	if (lowercase(key2)=="true") | (lowercase(key2)=="t")| (key2=="1")
		p = true
	elseif (lowercase(key2)=="false") | (lowercase(key2)=="f")| (key2=="0")
		p = false
	else
		return false
	end
	return key1==p
end
function comparekeys(key1::Number,key2::AbstractString)
	return key1==Meta.parse(key2)
end

function main(args)

    settings = ArgParseSettings(prog = "FITSexplore",
						 #version = @project_version,
						 version = "0.1",
						 add_version = true)

	settings.description =  "Simple tool to explore the content of FITS files.\n\n"*
							"Without any argument, it will display the name and the type of all HDU contained in the files TARGET."
    @add_arg_table! settings begin
		"--header", "-d"
			help = "header"
			action = :store_true
			help = "Print the whole FITS header."
        "--keyword", "-k"
			nargs = 1
			action = :append_arg
			arg_type = String
            help = "Print the value of the FITS header KEYWORD. This argument can be set multiple times to display several FITS keyword"
        "--filter", "-f"
            help = "filter"
			arg_type = String
			nargs = 2
			metavar = ["KEYWORD", "VALUE"]
			help = "Print all files where the FITS header KEYWORD = VALUE."
		"--recursive", "-r"
			help = "Recursively explore entire directories."
			action = :store_true
		"TARGET"
			nargs = '*'
			arg_type = String
            help = "List of all TARGET to explore. In conjunction with -r TARGET can contain directories."
			default = ["."]
    end

    parsed_args = parse_args(args, settings)

	args =parsed_args["TARGET"];

	files = Vector{String}()
	for arg in args
		if isdir(arg) && parsed_args["recursive"]
			files =  vcat(files,[root*"/"*filename for (root, dirs, TARGET) in walkdir(arg) for filename in TARGET  ])
		else
			files =  vcat(files,arg)
		end
	end

	head::Bool =  parsed_args["header"];

	if !isempty(parsed_args["keyword"])
		parse_keywords(files,parsed_args["keyword"])

	elseif !isempty(parsed_args["filter"])
		parse_filter(files,parsed_args["filter"])
	else
		for filename in files
			if isfile(filename)
				if endswith(filename,suffixes)
					if head
						@show read(FitsHeader,filename)
					else
						@show FITS(filename)
					end
				end
			end
		end
	end

end

end


text:
include("header.jl")

@testset "rnn" begin
    arrtype = KnetLayers.arrtype
    x    = arrtype(zeros(10,1))
    ind  = rand(1:10)
    x[ind,1] = 1.0
    l = LSTM(input=10,hidden=5,embed=5)
    l([ind]).y
    l(x).y
    l = SRNN(input=10,hidden=5,embed=5)
    l([ind]).y
    l(x).y
    l = GRU(input=10,hidden=5,embed=5)
    l([ind]).y
    l(x).y
    @test true
end


text:
module FAM

export coincidence, fa, quadgramscore, trigramscore, bigramscore

# include("Utils.jl")

using Main.UtilsM: removeExtras

# Index of Coincidence
function coincidence(text::String, alphabet::String)

    refined = removeExtras(text, alphabet)
    return sum([
        count(string(letter), refined) * (count(string(letter), refined) - 1)
        for letter in alphabet
    ]) / length(refined) / (length(refined) - 1) * length(alphabet)

end # function

# Basic Frequency Analysis
function fa(text::String, alphabet::String)

    refined = removeExtras(text, alphabet)
    # Count values in text.
    counted = Dict{Char, Float64}()
    for i in alphabet
        counted[i] = count(string(i), refined) / length(refined)
    end # for

    # Computer out:
    return counted

end # function

# Quadgram analysis - ONLY WORKS FOR ENGLISH
function quadgramScore(text::String)

    text = removeExtras(uppercase(text), "ABCDEFGHIJKLMNOPQRSTUVWXYZ")

    # Read the whole file into a dictionary for easy access
    qfdict = Dict()
    open("programs\\texts\\quadgrams.txt", "r") do qfreq

        while !eof(qfreq)

            temp = readline(qfreq)
            qfdict[temp[1:4]] = parse(Int, temp[6:end])

        end

    end

    # Score it
    return sum([
        get(qfdict, text[i-3:i], 0) for i in range(1, length = length(text)) if i > 3
    ])

end

# Trigram analysis - ONLY WORKS FOR ENGLISH
function trigramScore(text::String)

    text = removeExtras(uppercase(text), "ABCDEFGHIJKLMNOPQRSTUVWXYZ")

    # Read the whole file into a dictionary for easy access
    qfdict = Dict()
    open("programs\\texts\\trigrams.txt", "r") do qfreq

        while !eof(qfreq)

            temp = readline(qfreq)
            qfdict[temp[1:3]] = parse(Int, temp[5:end])

        end

    end

    # Score it
    return sum([
        get(qfdict, text[i-2:i], 0) for i in range(1, length = length(text)) if i > 2
    ])

end

# Bigram analysis - ONLY WORKS FOR ENGLISH
function bigramScore(text::String)

    text = removeExtras(uppercase(text), "ABCDEFGHIJKLMNOPQRSTUVWXYZ")

    # Read the whole file into a dictionary for easy access
    qfdict = Dict()
    open("programs\\texts\\bigrams.txt", "r") do qfreq

        while !eof(qfreq)

            temp = readline(qfreq)
            qfdict[temp[1:2]] = parse(Int, temp[4:end])

        end

    end

    # Score it
    return sum([
        get(qfdict, text[i-1:i], 0) for i in range(1, length = length(text)) if i > 1
    ])

end

end # module

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text:
function rowweights(array::Array{T})::Array{T} where T <: Number
    [sum(array[1:2:end]), sum(array[2:2:end])]
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using CountClick

push!(ARGS, "arg")
CountClick.julia_main()

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
"""
A simple simulation environment for a service mesh
"""
mutable struct ServiceMeshEnv{AT, R<:AbstractRNG} <: AbstractEnv
    # Mesh
    microservices::Vector{Microservice}

    # Jobs
    jobtypes::Vector{JobParams}

    # Track job data
    arrivals::Vector{Int}
    finished_jobs::Vector{Int}
    missed_deadlines::Vector{Int}
    dropped_jobs::Vector{Int}
    processed_jobs::Vector{Int}
    booting::Int
    closing::Int

    # Time
    time::Float64
    steplength::Float64 # Will be rounded down to even number of dts approximately
    dt::Float64

    # Spaces
    action_space::Space{Vector{AT}}
    state_space::Space

    # Other
    rng::R
end

function ServiceMeshEnv(;
        microservices::Int,
        min_scale::Int=1,
        max_scale::Int,
        boot_time::Float64,
        close_time::Float64 = boot_time,
        max_queue::Int,
        instance_cost::Float64,
        jobtypes::Vector{JobParams},
        steplength::Float64 = 1.0,
        dt::Float64 = 0.01,
        seed = 37,
        relative_action::Int=0,
        kwargs...
    )
    action_space = if relative_action != 0
        Space([-relative_action:relative_action for i in 1:microservices])
    else
        Space(ClosedInterval[min_scale..max_scale for i in 1:microservices])
    end
    ServiceMeshEnv(
        # Mesh
        [Microservice(
            boot_time = boot_time,
            close_time = close_time,
            min_scale = min_scale,
            max_scale = max_scale,
            max_queue = max_queue,
            instance_cost = instance_cost,
        ) for i in 1:microservices],

        # Job
        jobtypes,
        zeros(Int, length(jobtypes)),

        # Track job data
        zeros(Int, length(jobtypes)),
        zeros(Int, length(jobtypes)),
        zeros(Int, length(jobtypes)),
        zeros(Int, microservices),
        0,
        0,

        # Time
        0.0,
        steplength,
        dt,
        action_space,
        # Space(ClosedInterval[-1..1 for i in 1:microservices]),
        Space(
            ClosedInterval{Float64}[
                [min_scale..max_scale for i in 1:microservices]; # Running
                [0..1 for i in 1:microservices]; # Utilization
                # [0..max_queue for i in 1:microservices]; # Queue
                # [1..max_scale for i in 1:microservices]; # Target
                # [0..max_scale for i in 1:microservices]; # Booting
                # [0..max_scale for i in 1:microservices]; # Closing
                # [interval(jobtype.arrival) for jobtype in jobtypes]; # Arrival
            ]
        ),
        StableRNG(seed),
    )
end

"""
When a job finished execution on one node this finds the next step for it (next node, dropped, done)
"""
function propagate!(env::ServiceMeshEnv, job::Job)
    job.pathindex += 1
    if job.pathindex <= length(env.jobtypes[job.typeindex].path)
        dest = env.jobtypes[job.typeindex].path[job.pathindex]
        if length(env.microservices[dest].queue) < env.microservices[dest].running_nodes + env.microservices[dest].max_queue
            push!(env.microservices[dest].queue, job)
            job.countdown = env.jobtypes[job.typeindex].time[job.pathindex]
        else
            env.dropped_jobs[job.typeindex] += 1
        end
    else
        if env.time - job.arrivaltime < env.jobtypes[job.typeindex].deadline
            env.finished_jobs[job.typeindex] += 1
        else
            env.missed_deadlines[job.typeindex] += 1
        end
    end
end

Random.seed!(env::ServiceMeshEnv, seed) = Random.seed!(env.rng, seed)
RLBase.is_terminated(env::ServiceMeshEnv) = false

function RLBase.reset!(env::ServiceMeshEnv)
    for ms in env.microservices
        ms.target_nodes = 1
        ms.running_nodes = 1
        empty!(ms.booting_nodes)
        empty!(ms.closing_nodes)
        empty!(ms.queue)
    end

    fill!(env.arrivals, 0)
    fill!(env.finished_jobs, 0)
    fill!(env.missed_deadlines, 0)
    fill!(env.dropped_jobs, 0)
    fill!(env.processed_jobs, 0)
    env.booting = 0
    env.closing = 0

    env.time = 0
end

RLBase.action_space(env::ServiceMeshEnv) = env.action_space
RLBase.state_space(env::ServiceMeshEnv) = env.state_space

function RLBase.state(env::ServiceMeshEnv)
    s = vcat(
        map(ms -> ms.running_nodes, env.microservices),
        map(ms -> ms.running_nodes == 0 ? (isempty(ms.queue) ? 1.0 : ms.max_queue / 1)  : length(ms.queue) / ms.running_nodes, env.microservices),
        # begin
        #     tmp = map(ms->length(ms.queue), env.microservices)
        #     [tmp[1]; max.(tmp[1:end-1], tmp[2:end])]
        # end,
        # map(ms->length(ms.queue), env.microservices),
        # map(ms->ms.target_nodes, env.microservices),
        # map(ms->length(ms.booting_nodes), env.microservices),
        # map(ms->length(ms.closing_nodes), env.microservices),
        # env.arrivals,
    )

    return s
end

function RLBase.reward(env::ServiceMeshEnv)
    value = sum((env.finished_jobs[i] + env.jobtypes[i].missed_value_fraction * env.missed_deadlines[i]) * env.jobtypes[i].value for i in eachindex(env.finished_jobs))
    cost = sum((ms.running_nodes + length(ms.booting_nodes) + length(ms.closing_nodes)) * ms.instance_cost for ms in env.microservices)
    return value - cost
end

function (env::ServiceMeshEnv{<:ClosedInterval})(a::Vector{Int})
    @assert a in action_space(env)
    _step!(env, a)
end

# For relative actions
function (env::ServiceMeshEnv{<:UnitRange})(a::Vector{Int})
    for (i, ms) in enumerate(env.microservices)
        a[i] = clamp(ms.target_nodes + a[i], ms.min_scale, ms.max_scale)
    end

    @assert all(minlim .<= a .<= maxlim)
    _step!(env, a)
end

function _step!(env::ServiceMeshEnv, a)
    # Reset all counting
    fill!(env.finished_jobs, 0)
    fill!(env.missed_deadlines, 0)
    fill!(env.dropped_jobs, 0)
    fill!(env.processed_jobs, 0)
    fill!(env.arrivals, 0)
    env.booting = 0
    env.closing = 0

    # Update target nodes and boot/close towards this, should only be needed before loop
    scale.(env.microservices, a)

    for _ in 1:round(Int, env.steplength / env.dt)
        done_jobs = Job[]
        # Simulate each MS
        for (i, ms) in enumerate(env.microservices)
            # Work on jobs
            done_idx = Int[]
            for i in 1:min(ms.running_nodes, length(ms.queue))
                ms.queue[i].countdown -= env.dt
                if ms.queue[i].countdown <= 0
                    push!(done_idx, i)
                end
            end
            append!(done_jobs, ms.queue[done_idx])
            deleteat!(ms.queue, done_idx)
            env.processed_jobs[i] += length(done_idx)

            # Log here to catch data before done
            env.booting += length(ms.booting_nodes)
            env.closing += length(ms.closing_nodes)

            # Booting
            ms.booting_nodes .-= env.dt
            nb = length(ms.booting_nodes)
            filter!(>(0), ms.booting_nodes)
            ms.running_nodes += nb - length(ms.booting_nodes)
            # Closing
            ms.closing_nodes .-= env.dt
            filter!(>(0), ms.closing_nodes)
        end
        for job in done_jobs
            propagate!(env, job)
        end

        # These are the checked types of jobs when arriving which we can use in state
        new_arrivals = round.(Int, [get_arrivals(env.rng, jobtype.arrival, env.time) for jobtype in env.jobtypes])
        env.arrivals .+= new_arrivals

        # Adds arrivals to queue and shuffle before propagation
        newjobs = vcat((fill(i, n) for (i, n) in enumerate(new_arrivals))...)
        for jobtypeid in shuffle!(env.rng, newjobs)
            job = Job(jobtypeid, 0, 0, env.time)
            propagate!(env, job)
        end

        env.time += env.dt
    end
end

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
module BiMaps

export BiMap

importall Base

include("internal.jl")

type BiMap{L,R}
    left::BiMapInternal{true,L,R}
    right::BiMapInternal{false,L,R}
    function BiMap()
        l = Dict{L,R}()
        r = Dict{R,L}()
        left = BiMapInternal{true,L,R}(l, r)
        right = BiMapInternal{false,L,R}(l, r)
        @assert left.left === right.left
        @assert left.right === right.right
        new(left, right)
    end
end

function BiMap{L,R}(pairs::AbstractArray{(L,R)})
    bimap = BiMap{L,R}()
    for (lkey, rkey) in pairs
        bimap.left[lkey] = rkey
    end
    bimap
end

length(bimap::BiMap) = length(bimap.left)
isempty(bimap::BiMap) = isempty(bimap.left)

function show{L,R}(io::IO, bimap::BiMap{L,R})
    n = length(bimap)
    println(io, "BiMap{$L,$R} with $n ", n == 1 ? "entry" : "entries")
    for lkey in keys(bimap.left.left)
        rkey = bimap.left[lkey]
        write(io, "  ")
        show(io, lkey)
        write(io, " <=> ")
        show(io, rkey)
        write(io, '\n')
    end
end

start(bimap::BiMap) = start(bimap.left)
done(bimap::BiMap, s) = done(bimap.left, s)
next(bimap::BiMap, s) = next(bimap.left, s)

end # module