text:
# common args for all dmr tests
# to get exactly Rdistrom's run use λminratio=0.01 because our gamlr's have different defaults
testargs = Dict(:γ=>γdistrom, :λminratio=>0.01, :verbose=>false,:showwarnings=>true)

@testset "dmr" begin

f = @model(c ~ x + z + cat + y)
@test_show f "1-part model: [c ~ x + z + cat + y]"

dmrcoefs = dmr(covars, counts; testargs...)
coefs = coef(dmrcoefs)
@test size(coefs) == (p+1, d)
@test_throws MethodError coef(dmrcoefs, AllSeg())
@test coefs ≈ coefsRdistrom rtol=rtol
# println("rdist(coefs,coefsRdistrom)=$(rdist(coefs,coefsRdistrom))")

# test Int matrix for counts
dmrcoefsint = dmr(covars, countsint; testargs...)
coefsint = coef(dmrcoefsint)
@test coefsint == coefs

# test standardize=false
dmrcoefsunstd = dmr(covars, counts; standardize=false, testargs...)
coefsunstd = coef(dmrcoefsunstd)
@test coefsunstd != coefs
dmrcoefsunstd2 = dmr(covars, counts; standardize=false, local_cluster=false, testargs...)
coefsunstd2 = coef(dmrcoefsunstd2)
@test coefsunstd ≈ coefsunstd2

dmrcoefsb = fit(DMRCoefs, covars, counts; testargs...)
@test coef(dmrcoefsb) == coefs
dmrb = fit(DMR, covars, counts; testargs...)
@test coef(dmrb) == coefs
@test n == nobs(dmrcoefsb)
@test d == ncategories(dmrcoefsb)
@test p == ncovars(dmrcoefsb)

# select=MinBIC()
dmrcoefsb = fit(DMRCoefs, covars, counts; select=MinBIC(), testargs...)
@test coef(dmrcoefsb) != coefs
dmrb = fit(DMR, covars, counts; select=MinBIC(), testargs...)
@test coef(dmrb) != coefs

# serial run
dmrcoefss = dmr(covars, counts; parallel=false, testargs...)
@test coefs ≈ coef(dmrcoefss)

# using a dataframe and list of variables
dmrcoefsdf = fit(DMRCoefs, f, covarsdf, counts; testargs...)
@test coef(dmrcoefsdf) == coefs
dmrpathsdf = fit(DMRPaths, f, covarsdf, counts; testargs...)
@test coef(dmrpathsdf) ≈ coefs

dmrcoefs2 = dmr(covars, counts; local_cluster=false, testargs...)
coefs2 = coef(dmrcoefs2)
@test coefs ≈ coefs2
dmrcoefs2 = fit(DMRCoefs, covars, counts; local_cluster=false, testargs...)
@test coef(dmrcoefs2) == coefs2

# keeping regulatrization paths
dmrPaths = dmrpaths(covars, counts; testargs...)
dmrPaths2 = fit(DMRPaths, covars, counts; testargs...)
# serial paths
dmrPathss = fit(DMRPaths, covars, counts; parallel=false, testargs...)
@test coef(dmrPaths) == coef(dmrPaths2)
@test coef(dmrPaths) ≈ coefs
@test coef(dmrPathss) ≈ coefs
coefsall = coef(dmrPaths, AllSeg())
@test size(coefsall,1) > 1
@test size(coefsall,2) == p+1
@test size(coefsall,3) == d
@test n == nobs(dmrPaths)
@test d == ncategories(dmrPaths)
@test p == ncovars(dmrPaths)

z = srproj(coefs, counts)
zb = srproj(dmrcoefs, counts)
zc = srproj(dmrPaths, counts)
@test z == zb
@test z ≈ zc
@test size(z) == (size(counts,1),p+1)

regdata = DataFrame(y=covars[:,1], z=z[:,1], m=z[:,2])
lm1 = lm(@formula(y ~ z+m), regdata)
r21 = adjr2(lm1)

η = predict(dmrPaths,newcovars)
@test η ≈ predictRdistrom rtol=5rtol
# rdist(η,predictRdistrom)
@test sum(η, dims=2) ≈ ones(size(newcovars, 1))

@test_throws ErrorException predict(dmrcoefs,newcovars)

@test z ≈ zRdistrom rtol=rtol
# println("rdist(z,zRdistrom)=$(rdist(z,zRdistrom))")

regdata3 = DataFrame(y=covars[:,1], z=zRdistrom[:,1], m=zRdistrom[:,2])
lm3 = lm(@formula(y ~ z+m), regdata3)
r23 = adjr2(lm3)

@test r21 ≈ r23 rtol=rtol
# println("rdist(r21,r23)=$(rdist(r21,r23))")

# project in a single direction
z1 = srproj(coefs, counts, projdir)
z1b = srproj(dmrcoefs, counts, projdir)
z1c = srproj(coefs, counts, projdir)
@test z1 ≈ z[:,[projdir,end]]

z1dense = srproj(coefs, Matrix(counts), projdir)
z1denseb = srproj(dmrcoefs, Matrix(counts), projdir)
z1densec = srproj(dmrPaths, Matrix(counts), projdir)
@test z1dense == z1denseb
@test z1dense ≈ z1densec
@test z1dense ≈ z1

@test z1 ≈ z1Rdistrom rtol=rtol
X1, X1_nocounts, inz = srprojX(coefs,counts,covars,projdir; includem=true)
@test X1_nocounts == [ones(n,1) covars[:,1:4]]
@test X1 == [X1_nocounts z1]
@test inz == [1]

X1b, X1_nocountsb, inz = srprojX(dmrcoefs,counts,covars,projdir; includem=true)
@test X1 == X1b
@test X1_nocounts == X1_nocountsb
@test inz == [1]

X2, X2_nocounts, inz = srprojX(coefs,counts,covars,projdir; includem=false)
@test X2_nocounts == X1_nocounts
@test X2 == X1[:,1:end-1]
@test inz == [1]

X2b, X2_nocountsb, inz = srprojX(dmrcoefs,counts,covars,projdir; includem=false)
@test X2 == X2b
@test X2_nocounts == X2_nocountsb
@test inz == [1]

# project in a single direction, focusing only on cj
focusj = 3
z1j = srproj(coefs, counts, projdir; focusj=focusj)
z1jdense = srproj(coefs, Matrix(counts), projdir; focusj=focusj)
@test z1jdense == z1j

X1j, X1j_nocounts, inz = srprojX(coefs,counts,covars,projdir; includem=true, focusj=focusj)
@test X1j_nocounts == [ones(n,1) covars[:,1:4]]
@test X1j == [X1_nocounts z1j]
@test inz == [1]

X2j, X2j_nocounts, inz = srprojX(coefs,counts,covars,projdir; includem=false, focusj=focusj)
@test X2j_nocounts == X1_nocounts
@test X2j == X1j[:,1:end-1]
@test inz == [1]

X1jfull, X1jfull_nocounts, inz = srprojX(coefs,Matrix(counts),covars,projdir; includem=true, focusj=focusj)
@test X1jfull ≈ X1j
@test X1jfull_nocounts ≈ X1jfull_nocounts
@test inz == [1]

X2jfull, X2jfull_nocounts, inz = srprojX(coefs,Matrix(counts),covars,projdir; includem=false, focusj=focusj)
@test X2jfull ≈ X2j
@test X2jfull_nocounts ≈ X2jfull_nocounts
@test inz == [1]

# MNIR
mnir = fit(CIR{DMR,LinearModel},covars,counts,projdir; nocounts=true, testargs...)
@test coefbwd(mnir) ≈ coef(dmrcoefs)
mnirglm = fit(CIR{DMR,GeneralizedLinearModel},covars,counts,projdir,Gamma(); nocounts=false, testargs...)
@test coefbwd(mnirglm) ≈ coef(dmrcoefs)
@test !(coeffwd(mnirglm) ≈ coeffwd(mnir))
@test_throws ErrorException coeffwd(mnirglm; nocounts=true)
@test !(predict(mnir,covars[1:10,:],counts[1:10,:]) ≈ predict(mnirglm,covars[1:10,:],counts[1:10,:]))
@test !(predict(mnir,covars[1:10,:],counts[1:10,:]) ≈ predict(mnir,covars[1:10,:],counts[1:10,:];nocounts=true))
@test_throws ErrorException predict(mnirglm,covars[1:10,:],counts[1:10,:];nocounts=true)

mnirdf = fit(CIR{DMR,LinearModel},f,covarsdf,counts,"y"; nocounts=true, testargs...)
@test coefbwd(mnirdf) ≈ coef(dmrcoefs)
@test coeffwd(mnirdf) ≈ coeffwd(mnir)
@test coeffwd(mnirdf) != coeffwd(mnir; nocounts=true)
@test coef(mnirdf) ≈ coef(mnir)
@test_show mnirdf string(f)
@test_show mnirdf "Forward model coefficients for predicting y"

mnirglmdf = fit(CIR{DMR,GeneralizedLinearModel},f,covarsdf,counts,:y,Gamma(); nocounts=false, testargs...)
@test coefbwd(mnirglmdf) ≈ coef(dmrcoefs)
@test_throws ErrorException coeffwd(mnirglmdf; nocounts=true)
@test_show mnirglmdf string(f)
@test_show mnirglmdf "Forward model coefficients for predicting y"

# select=MinBIC()
mnirdfb = fit(CIR{DMR,LinearModel},f,covarsdf,counts,:y; nocounts=true, select=MinBIC(), testargs...)
@test coefbwd(mnirdf) != coefbwd(mnirdfb)
@test coeffwd(mnirdf) != coeffwd(mnirdfb)

# #### debug start
# mm = mnirdf
# mdf = covarsdf[1:10,:]
# newTerms = StatsModels.dropresponse!(mm.mf.terms)
# # create new model frame/matrix
# newTerms.intercept = true
# mf = ModelFrame(newTerms, mdf; contrasts = mm.mf.contrasts)
# mf.terms.intercept = false
# newX = ModelMatrix(mf).m
# if !all(mf.nonmissing)
# counts = counts[mf.nonmissing,:]
# end
# yp = predict(mm, newX, counts; kwargs...)
# out = missings(eltype(yp), size(df, 1))
# out[mf.nonmissing] = yp
#
# predict(mnirdf,covarsdf[1:10,:],counts[1:10,:])
# size(covarsdf)
# mnirdf.model.projdir
# #### debug end
@test !(predict(mnirdf,covarsdf[1:10,:],counts[1:10,:]) ≈ predict(mnirglmdf,covarsdf[1:10,:],counts[1:10,:]))
@test !(predict(mnirdf,covarsdf[1:10,:],counts[1:10,:]) ≈ predict(mnirglmdf,covarsdf[1:10,:],counts[1:10,:]))
@test !(predict(mnirdf,covarsdf[1:10,:],counts[1:10,:]) ≈ predict(mnirdf,covarsdf[1:10,:],counts[1:10,:];nocounts=true))
@test_throws ErrorException predict(mnirglmdf,covarsdf[1:10,:],counts[1:10,:];nocounts=true)

zlm = lm(hcat(ones(n,1),z1,covars[:,1:4]),covars[:,projdir])
@test r2(zlm) ≈ r2(mnir)
@test adjr2(zlm) ≈ adjr2(mnir)
@test predict(zlm,hcat(ones(10,1),z1[1:10,:],covars[1:10,1:4])) ≈ predict(mnir,covars[1:10,:],counts[1:10,:])
@test predict(zlm,hcat(ones(10,1),z1[1:10,:],covars[1:10,1:4])) ≈ predict(mnirdf,covarsdf[1:10,:],counts[1:10,:])

zlmnocounts = lm(hcat(ones(n,1),covars[:,1:4]),covars[:,projdir])
@test r2(zlmnocounts) ≈ r2(mnir; nocounts=true)
@test adjr2(zlmnocounts) ≈ adjr2(mnir; nocounts=true)
@test predict(zlmnocounts,hcat(ones(10,1),covars[1:10,1:4])) ≈ predict(mnir,covars[1:10,:],counts[1:10,:]; nocounts=true)
@test predict(zlmnocounts,hcat(ones(10,1),covars[1:10,1:4])) ≈ predict(mnirdf,covarsdf[1:10,:],counts[1:10,:]; nocounts=true)

end

#########################################################################3
# degenerate cases
#########################################################################3

@testset "dmr degenerate cases" begin

@info("Testing dmr degenerate cases. The 2 following warnings by workers are expected ...")

f = @model(c ~ x + z + cat + y)

# always one (zero var) counts columns
zcounts = deepcopy(counts)
zcounts[:,2] = zeros(n)
zcounts[:,3] = ones(n)

# make sure we are not adding all zero obseravtions
m = sum(zcounts,dims=2)
@test sum(m .== 0) == 0

# this one should warn on dimension 2 but @test_warn doen't capture the workers' warnings
dmrzcoefs = dmr(covars, zcounts; testargs...)
zcoefs = coef(dmrzcoefs)
@test size(zcoefs) == (p+1, d)
@test zcoefs[:,2] ≈ zeros(p+1)

dmrzcoefs2 = dmr(covars, zcounts; local_cluster=false, testargs...)
zcoefs2 = coef(dmrzcoefs2)
@test zcoefs2 ≈ zcoefs

# serial runs can also test for warnings
dmrzcoefs3 = @test_logs (:warn, r"failed on count dimension 2") dmr(covars, zcounts; parallel=false, testargs...)
zcoefs2 = coef(dmrzcoefs2)
@test zcoefs2 ≈ zcoefs

dmrzpaths3 = @test_logs (:warn, r"failed for countsj") dmrpaths(covars, zcounts; parallel=false, testargs...)
zcoefs3 = coef(dmrzpaths3)
@test zcoefs3 ≈ zcoefs rtol=rtol

# test an observation with all zeros
zcounts = deepcopy(counts)
zcounts[1,:] .= 0.0
m = sum(zcounts, dims=2)
@test sum(m .== 0) == 1
dmrzcoefs = @test_logs (:warn, r"omitting 1") dmr(covars, zcounts; testargs...)
dmrzcoefs2 = dmr(covars[2:end,:], counts[2:end,:]; testargs...)
zcoefs3 = coef(dmrzcoefs)
@test size(zcoefs3) == (p+1, d)
@test nobs(dmrzcoefs) == n-1
@test zcoefs3 == coef(dmrzcoefs2)

zcovarsdf = deepcopy(covarsdf)
allowmissing!(zcovarsdf, :y)
zcovarsdf[1,:y] = missing
dmrzcoefsdf = fit(DMR, f, zcovarsdf, counts; testargs...)
zcoefsdf = coef(dmrzcoefsdf)
@test zcoefsdf == zcoefs3

zmnirdf = fit(CIR{DMR,LinearModel},f,zcovarsdf,counts,:y; nocounts=true, testargs...)
@test_show zmnirdf string(f)
@test_show zmnirdf "Forward model coefficients for predicting y"

# a missing lhs variable shouldn't matter
zyhat = predict(zmnirdf, zcovarsdf, counts)
@test !any(ismissing,zyhat)

# but a missing rhs variable means an observation is dropped
allowmissing!(zcovarsdf, :x)
zcovarsdf[1,:x] = missing
zyhat = predict(zmnirdf, zcovarsdf, counts)
@test ismissing(zyhat[1])
@test !any(ismissing,zyhat[2:end])

end


text:
# File for getting the data for Figure 4

using Pkg
using Statistics
using Random, Distributions
using SparseGrids
using BenchmarkTools
using DataFrames
using JLD2

include("/home/christian/UniMA/EMMA/src/RFR.jl")
include("/home/christian/UniMA/EMMA/src/cross_val.jl")
include("/home/christian/UniMA/EMMA/src/aux_functions.jl")

import Cairo, Fontconfig

# Set seed for reproducability
Random.seed!(68151)

# Set parameters
n = 2000
d = 10
σ = 1
m_features = d

# Make dictionaries for CART and weighted approach
d1 = Dict{Symbol, Vector{Float64}}(
    :max_features => [m_features],
    :n_trees => [30])

d2 = copy(d1)
d2[:α] = [0] # CART has α=0

# List of observations to go through
n_list = [250, 500, 1000, 2000, 4000, 8000, 16000]
res_mat = zeros(length(n_list), 2)

for (ind,n1) in enumerate(n_list)
    x_train, x_test, y_train, y_test = make_data(n1, d, "friedman", σ)
    cv = cross_val(d1, random_state=0)
    fit!(cv, x_train, y_train)

    # Get optimal parameter for α
    rf_α = best_model(cv, "mse")
    rf0 = RFR(param_dict = d2)

    fit!(rf_α, x_test, y_test)
    fit!(rf0, x_test, y_test)

    # Save average depth of all trees in forest
    res_mat[ind, 1] = average_depth(rf0)
    res_mat[ind, 2] = average_depth(rf_α)
    println(ind)
end

# Make dataframe and save data
df_tmp = DataFrame(res_mat, :auto)
rename!(df_tmp, ["unweighted", "weighted"])
df_tmp[!, :n] = n_list

jldsave("data/tree_depth.jld2"; df_tmp)

ext: jl

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
import Base.findin
import Base.insert!
import Base.push!

"""
Part of the general awesomeness of MathGL for scientific graphics is that it never does anything that it is not explicitly told to.
This allows for extreme flexibility in how your graphics look, but it also tends to create code that looks like this:
dat = generate_some_data()
gr = mglGraph()
SetRange(gr, 'x', size(dat, 1))
SetRange(gr, 'y', size(dat, 2))
SetRange(gr, 'z', minimum(dat), maximum(dat))
Surf(gr, data)
Axis(gr)
Box(gr)

rather than like this:

dat = generate_some_data()
Surf(gr, dat)

This is all very well, not to say entirely appropriate, when you have a program that generates lots of figures. It provides the ultimate control.
But it is a very different story at the REPL. To get a pretty picture, you have to type seven commands; If you mis-type number seven, then you have
to start all over again. But condensing these seven operations into one simple command, e.g. 'Surf', would create the opposite problem. You would
have to have a ton of different definitions of 'Surf' for however you wanted your graph to look.

The plotOpStack idea brings a new approach to this situation. There are two different methods for Surf: One follows the MathGL API faithfully to enable full
control by you when you need it. The other returns a set of operations to be performed on an mglGraph, stored in an array along with descriptive names and boolean
switches to turn them on and off. Thus the user can type:
#At the beginning of the session, e.g. in the user's julia startup file
import Base.show
function show(io::IO, ops::plotOpStack)
	ShowImage(draw(ops), "FancyImageViewer"))
end

# At the REPL
dat = generate_some_data()
ops = Surf(dat) #Note: when the first parameter is not an mglGraph object, the method that returns a plotOpStack is called

Then the plot will be shown. If the user wishes to change some operation, he/she can easily disable things by setting elements of ops.switches to false
or call insert!(ops, someDrawingFunction, someIndex) too add the operation someDrawingFunction to the stack at position someIndex
"""
struct plotOpStack
    ops::Array{Function, 1}
    names::Array{String, 1}
    switches::Array{Bool, 1}
    plotOpStack() = new(Array{Function, 1}(0), Array{String, 1}(0), Array{Bool, 1}(0))
end

"""Push operation _op_ with name _name_ and enabled state _on_ onto the end of a plotOpStack"""
function push!(opStack::plotOpStack, op::Function, name::String="", on::Bool=true)
	push!(opStack.ops, op)
	push!(opStack.names, name)
	push!(opStack.switches, on)
end

"""Insert operation _op_ into the plotOpStack at index _index_ with name _name_ and enabled/disabled state given by _on_"""
function insert!(opStack::plotOpStack, index::Integer, op::Function, name::String="", on::Bool=true)
	insert!(opStack.ops, index, op)
	insert!(opStack.names, index, name)
	insert!(opStack.switches, index, on)
end

"""Return the index(es) of the operation(s) with name(s) `name` in the plotOpStack `opStack`"""
function findin(opStack::plotOpStack, name::String)
    return find(x->x==name, opStack.names)
end

"""Disable the operation at index(es) `indecies...` in the plotOpStack `opStack`"""
function disable!(opStack::plotOpStack, indecies::Integer...)
    for index in indecies
	opStack.switches[index] = false
    end
    return nothing
end

"""Disable the operation named `name` in the plotOpStack `opStack`, if it exists"""
function disable!(opStack::plotOpStack, name::String)
    indecies = findin(opStack, name)
    if length(indecies) == 0
	warn("Cannot disable operation $name because it is not present in the plotOpStack")
    end
    disable!(opStack, indecies...)
end

"""Apply the functions in _opStack_ to _gr_ in order"""
function draw!(gr::mglGraph, opStack::plotOpStack)
    nOps = size(opStack.ops, 1)
    for i in 1:nOps
	if opStack.switches[i]
		opStack.ops[i](gr)
		if isinteractive()
#			println("Drawing ", opStack.names[i])
		end
	end
    end
end

"""Create a mglGraph and apply the operations in _opStack_ to it"""
function draw(opStack::plotOpStack, width::Integer=-1, height::Integer=-1)
	if width == -1
		width = default_width
	end
	if height == -1
		height = default_height
	end

	gr = mglGraph(width, height)
	draw!(gr, opStack)
	return gr
end

export push!
export findin
export insert!
export disable!

export plotOpStack


text:
module Resources

include("meta.jl")
export TypeMeta, ObjectMeta

end # module

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
using Test

# Build the example docs
include("examples/make.jl")

# Test missing docs
include("missingdocs/make.jl")

# Primary @testset

# Error reporting.
println("="^50)
@info("The following errors are expected output.")
include(joinpath("errors", "make.jl"))
@info("END of expected error output.")
println("="^50)

@testset "Documenter" begin
    # Unit tests for module internals.
    include("utilities.jl")

    # NavNode tests.
    include("navnode.jl")

    # DocSystem unit tests.
    include("docsystem.jl")

    # DocTest unit tests.
    include("doctests/doctests.jl")

    # DOM Tests.
    include("dom.jl")

    # MDFlatten tests.
    include("mdflatten.jl")

    # HTMLWriter
    include("htmlwriter.jl")

    # Mock package docs.
    include("examples/tests.jl")

    # Documenter package docs with other formats.
    include("formats/markdown.jl")
    include("formats/latex.jl")

    # A simple build outside of a Git repository
    include("nongit/tests.jl")
end

# Additional tests

## `Markdown.MD` to `DOM.Node` conversion tests.
module MarkdownToNode
    import Documenter.DocSystem
    import Documenter.Writers.HTMLWriter: mdconvert

    # Exhaustive Conversion from Markdown to Nodes.
    for mod in Base.Docs.modules
        for (binding, multidoc) in DocSystem.getmeta(mod)
            for (typesig, docstr) in multidoc.docs
                md = DocSystem.parsedoc(docstr)
                string(mdconvert(md))
            end
        end
    end
end


text:
__precompile__(true)

module Tri_Tet_Intersections

using StaticArrays
using ForwardDiff: value
using LinearAlgebra
using NumericalTricks
using RigidBodyDynamics.Spatial

include("poly_eight.jl")
include("utility.jl")
include("quadrature.jl")
include("static_clip.jl")
include("plane_tet_intersection.jl")
include("test_utility.jl")

export
    # poly_eight.jl
    poly_eight,
    zero_small_coordinates,

    # utility.jl
    asMatOnePad,
    asMat,

    # quadrature.jl
    TriTetQuadRule,
    getTriQuadRule,
    getTetQuadRule,

    # static_clip.jl
    poly_eight,
    clip_in_tet_coordinates,

    # plane_tet_intersection.jl
    clip_plane_tet,

    # test_utility.jl
    roll_non_degenerate_tet,
    make_4_sided,
    dist_from_plane,
    project_into_plane

end


text:
abstract type AbstractScalarFunction <: JuMP.AbstractJuMPScalar end

struct Volume{S} <: AbstractScalarFunction
    set::S
end
Polyhedra.volume(variable::Union{SetVariableRef, Sets.Projection{<:SetVariableRef}, Sets.AbstractSet{<:JuMP.AbstractJuMPScalar}}) = Volume(variable)

function JuMP.set_objective(model::JuMP.Model,
                            sense::MOI.OptimizationSense,
                            func::AbstractScalarFunction)
    d = data(model)
    @assert d.state == Modeling
    d.objective_sense = sense
    d.objective = func
end

function load(model::JuMP.Model, f::AbstractScalarFunction)
    JuMP.set_objective_sense(model, objective_sense(model, f))
    JuMP.set_objective_function(model, objective_function(model, f))
end

include("affine.jl")
include("nth_root.jl")
include("L1_heuristic.jl")


text:
function helper_fn(z, attempt)
	single_data_name = z
	proc_id = myid();

	helper_start_time = time()

	train_transaction_counts = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_transaction_counts.csv") ) )
	train_assortments = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_assortments/", single_data_name, "_assortments.csv") ) )

	test_transaction_counts = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_R2_temporally_split/week3_to_week6/", single_data_name, "_transaction_counts.csv") ) )
	test_assortments = convert(Array{Int64}, readcsv( string("../tcm_data/IRI_data_R2_temporally_split/week3_to_week6/", single_data_name, "_assortments.csv") ) )

	num_train_transactions = sum(train_transaction_counts)
	num_test_transactions = sum(test_transaction_counts)

	num_train_assortments = size(train_assortments, 1)
	num_test_assortments = size(test_assortments, 1)

	N = 10

	time_limit = 1200

	perTol = 1e-6
	perTol_EM = 1e-5;
	time_limit_EM = time_limit;
	K_sample = 2000;

	## RANKING-BASED MODEL ##
	ranking_lambda, orderings, ranking_loglik, ranking_elapsed_time = tcm_estimateRanking(N, train_assortments, train_transaction_counts, time_limit, perTol)
	ranking_K = length(ranking_lambda)

	test_ranking_predict, test_ranking_predict_by_class = tcm_predictRanking(N, ranking_K, ranking_lambda, orderings, test_assortments)
	ranking_KL, ranking_noinf_count, ranking_KL_avg = tcm_evaluateKL(N, test_ranking_predict, test_transaction_counts)

	avg_CS_size, weighted_avg_CS_size, max_CS_size = tcm_rankingStats(N, ranking_lambda, orderings)

	row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"ranking", ",",
							ranking_loglik, ",",
							ranking_K, ",",
							perTol, ",",
							ranking_KL, ",",
							ranking_noinf_count, ",",
							ranking_KL_avg, ",",
							ranking_elapsed_time, ",",
							avg_CS_size, ",",
							weighted_avg_CS_size, ",",
							max_CS_size)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(ranking_outcsvhandle, row_string);
		flush(ranking_outcsvhandle);
	end

	## FOREST, RTS, COLD-START ##

	depth_list = [3,4,5,6,7]

	# initial_left, initial_right, initial_product, initial_isLeaf = tcm_createBasicForest(N)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_loglik, forest_elapsed_time = tcm_estimateForest_RTS(N, train_assortments, train_transaction_counts, K_sample, depth_limit, perTol_EM, time_limit_EM)
		K = length(lambda)

		test_forest_predict, test_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, test_assortments)
		forest_KL, forest_noinf_count, forest_KL_avg = tcm_evaluateKL(N, test_forest_predict, test_transaction_counts)

		simple_avg_leaves, weighted_avg_leaves, max_leaves = tcm_leafStats(lambda, forest_isLeaf)
		simple_avg_depth, weighted_avg_depth, max_depth = tcm_depthStats(lambda, forest_left, forest_right, forest_isLeaf)

		row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"forest-cold-RTS", ",",
							depth_limit, ",",
							K_sample, ",",
							perTol_EM, ",",
							forest_loglik, ",",
							K, ",",
							perTol, ",",
							forest_KL, ",",
							forest_noinf_count, ",",
							forest_KL_avg, ",",
							forest_elapsed_time, ",",
							simple_avg_leaves, ",",
							weighted_avg_leaves, ",",
							max_leaves, ",",
							simple_avg_depth, ",",
							weighted_avg_depth, ",",
							max_depth)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end
	end

	## FOREST, WARM-START ##
	depth_list = [3,4,5,6,7]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_convertRankingToForest(N, orderings)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_loglik, forest_elapsed_time = tcm_estimateForest_RTS_withInitial(N, train_assortments, train_transaction_counts, K_sample, depth_limit, perTol_EM, time_limit_EM,
																												initial_left, initial_right, initial_product, initial_isLeaf)
		K = length(lambda)

		test_forest_predict, test_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, test_assortments)
		forest_KL, forest_noinf_count, forest_KL_avg = tcm_evaluateKL(N, test_forest_predict, test_transaction_counts)

		simple_avg_leaves, weighted_avg_leaves, max_leaves = tcm_leafStats(lambda, forest_isLeaf)
		simple_avg_depth, weighted_avg_depth, max_depth = tcm_depthStats(lambda, forest_left, forest_right, forest_isLeaf)

		row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							num_train_transactions, ",",
							num_test_transactions, ",",
							proc_id, ",",
							attempt, ",",
							"forest-ranking-warm-RTS", ",",
							depth_limit, ",",
							K_sample, ",",
							perTol_EM, ",",
							forest_loglik, ",",
							K, ",",
							perTol, ",",
							forest_KL, ",",
							forest_noinf_count, ",",
							forest_KL_avg, ",",
							forest_elapsed_time, ",",
							simple_avg_leaves, ",",
							weighted_avg_leaves, ",",
							max_leaves, ",",
							simple_avg_depth, ",",
							weighted_avg_depth, ",",
							max_depth)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end
	end

	helper_end_time = time()

	return helper_end_time - helper_start_time;

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
"""
```
JSON.lower(
    t::T,
) where {
    T <: Union{
        AbstractCrystalStruct,
        AbstractMobility,
        AbstractIntegrator,
        AbstractDlnSeg,
        AbstractDlnStr,
        AbstractDistribution,
    },
}

JSON.lower(t::nodeTypeDln)
```
Extensions to `JSON.lower` for custom types. Allows these variables to be serialised properly.
"""
function JSON.lower(
    t::T,
) where {
    T <: Union{
        AbstractCrystalStruct,
        AbstractMobility,
        AbstractIntegrator,
        AbstractDlnSeg,
        AbstractDlnStr,
        AbstractDistribution,
        DispatchRegularCuboidMesh,
        LinearElement,
        CantileverLoad,
    },
}
    return string(t)
end
JSON.lower(t::nodeTypeDln) = Int(t)
JSON.lower(t::nodeTypeFE) = Int(t)
JSON.lower(t::Ptr{SuiteSparse.CHOLMOD.C_Factor{Float64}}) = nothing
"""
```
saveJSON(filename::AbstractString, args...; mode::AbstractString = "w")
```
Wrapper for `JSON.print`.
"""
function saveJSON(filename::AbstractString, args...; mode::AbstractString = "w")
    open(filename, mode) do io
        return length(args) == 1 ? JSON.print(io, args...) : JSON.print(io, args)
    end
    return nothing
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
module GLPKInterfaceLP

import GLPK
importall MathProgBase.SolverInterface
importall ..GLPKInterfaceBase

export
    GLPKSolverLP,
    optimize!,
    loadproblem!,
    writeproblem,
    getvarLB,
    setvarLB!,
    getvarUB,
    setvarUB!,
    getconstrLB,
    setconstrLB!,
    getconstrUB,
    setconstrUB!,
    getobj,
    setobj!,
    addvar!,
    addconstr!,
    setsense!,
    getsense,
    numvar,
    numconstr,
    status,
    getobjval,
    getsolution,
    getconstrsolution,
    getreducedcosts,
    getconstrduals,
    getrawsolver,
    getinfeasibilityray,
    getunboundedray

type GLPKMathProgModelLP <: GLPKMathProgModel
    inner::GLPK.Prob
    method::Symbol
    param::Union{GLPK.SimplexParam, GLPK.InteriorParam}
end

type GLPKSolverLP <: AbstractMathProgSolver
    presolve::Bool
    method::Symbol
    opts
    function GLPKSolverLP(;presolve::Bool=false, method::Symbol=:Simplex, opts...)
        method in [:Simplex, :Exact, :InteriorPoint] ||
            error("""
                  Unknown method for GLPK LP solver: $method
                         Allowed methods:
                           :Simplex
                           :Exact
                           :InteriorPoint""")
        new(presolve, method, opts)
    end
end

function LinearQuadraticModel(s::GLPKSolverLP)
    if s.method == :Simplex || s.method == :Exact
        param = GLPK.SimplexParam()
        if s.presolve
            param.presolve = GLPK.ON
        end
    elseif s.method == :InteriorPoint
        param = GLPK.InteriorParam()
        if s.presolve
            warn("Ignored option: presolve")
        end
    else
        error("This is a bug")
    end
    param.msg_lev = GLPK.MSG_ERR
    for (k,v) in s.opts
        i = findfirst(x->x==k, fieldnames(typeof(param)))
        if i > 0
            t = typeof(param).types[i]
            setfield!(param, i, convert(t, v))
        else
            warn("Ignored option: $(string(k))")
        end
    end
    lpm = GLPKMathProgModelLP(GLPK.Prob(), s.method, param)
    return lpm
end

function optimize!(lpm::GLPKMathProgModelLP)
    if lpm.method == :Simplex
        solve = GLPK.simplex
    elseif lpm.method == :Exact
        solve = GLPK.exact
    elseif lpm.method == :InteriorPoint
        solve = GLPK.interior
    else
        error("bug")
    end
    return solve(lpm.inner, lpm.param)
end

function status(lpm::GLPKMathProgModelLP)
    if lpm.method == :Simplex || lpm.method == :Exact
        get_status = GLPK.get_status
    elseif lpm.method == :InteriorPoint
        get_status = GLPK.ipt_status
    else
        error("bug")
    end
    s = get_status(lpm.inner)
    if s == GLPK.OPT
        return :Optimal
    elseif s == GLPK.INFEAS
        return :Infeasible
    elseif s == GLPK.UNBND
        return :Unbounded
    elseif s == GLPK.FEAS
        return :Feasible
    elseif s == GLPK.NOFEAS
        return :Infeasible
    elseif s == GLPK.UNDEF
        return :Undefined
    else
        error("Internal library error")
    end
end

function getobjval(lpm::GLPKMathProgModelLP)
    if lpm.method == :Simplex || lpm.method == :Exact
        get_obj_val = GLPK.get_obj_val
    elseif lpm.method == :InteriorPoint
        get_obj_val = GLPK.ipt_obj_val
    else
        error("bug")
    end
    return get_obj_val(lpm.inner)
end

function getsolution(lpm::GLPKMathProgModelLP)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)

    if lpm.method == :Simplex || lpm.method == :Exact
        get_col_prim = GLPK.get_col_prim
    elseif lpm.method == :InteriorPoint
        get_col_prim = GLPK.ipt_col_prim
    else
        error("bug")
    end

    x = Array(Float64, n)
    for c = 1:n
        x[c] = get_col_prim(lp, c)
    end
    return x
end

function getconstrsolution(lpm::GLPKMathProgModelLP)
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)

    if lpm.method == :Simplex || lpm.method == :Exact
        get_row_prim = GLPK.get_row_prim
    elseif lpm.method == :InteriorPoint
        get_row_prim = GLPK.ipt_row_prim
    else
        error("bug")
    end

    x = Array(Float64, m)
    for r = 1:m
        x[r] = get_row_prim(lp, r)
    end
    return x
end

function getreducedcosts(lpm::GLPKMathProgModelLP)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)

    if lpm.method == :Simplex || lpm.method == :Exact
        get_col_dual = GLPK.get_col_dual
    elseif lpm.method == :InteriorPoint
        get_col_dual = GLPK.ipt_col_dual
    else
        error("bug")
    end

    x = Array(Float64, n)
    for c = 1:n
        x[c] = get_col_dual(lp, c)
    end
    return x
end

function getconstrduals(lpm::GLPKMathProgModelLP)
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)

    if lpm.method == :Simplex || lpm.method == :Exact
        get_row_dual = GLPK.get_row_dual
    elseif lpm.method == :InteriorPoint
        get_row_dual = GLPK.ipt_row_dual
    else
        error("bug")
    end

    x = Array(Float64, m)
    for r = 1:m
        x[r] = get_row_dual(lp, r)
    end
    return x
end

# The functions getinfeasibilityray and getunboundedray are adapted from code
# taken from the LEMON C++ optimization library. This is the copyright notice:
#
### Copyright (C) 2003-2010
### Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
### (Egervary Research Group on Combinatorial Optimization, EGRES).
###
### Permission to use, modify and distribute this software is granted
### provided that this copyright notice appears in all copies. For
### precise terms see the accompanying LICENSE file.
###
### This software is provided "AS IS" with no warranty of any kind,
### express or implied, and with no claim as to its suitability for any
### purpose.

function getinfeasibilityray(lpm::GLPKMathProgModelLP)
    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
    elseif lpm.method == :InteriorPoint
        error("getinfeasibilityray is not available when using the InteriorPoint method")
    else
        error("bug")
    end

    m = GLPK.get_num_rows(lp)

    ray = zeros(m)

    ur = GLPK.get_unbnd_ray(lp)
    if ur != 0
        if ur <= m
            k = ur
            get_stat = GLPK.get_row_stat
            get_bind = GLPK.get_row_bind
            get_prim = GLPK.get_row_prim
            get_ub = GLPK.get_row_ub
        else
            k = ur - m
            get_stat = GLPK.get_col_stat
            get_bind = GLPK.get_col_bind
            get_prim = GLPK.get_col_prim
            get_ub = GLPK.get_col_ub
        end

        get_stat(lp, k) == GLPK.BS || error("unbounded ray is primal (use getunboundedray)")

        ray[get_bind(lp, k)] = (get_prim(lp, k) > get_ub(lp, k)) ? -1 : 1

        GLPK.btran(lp, ray)
    else
        eps = 1e-7
        for i = 1:m
            idx = GLPK.get_bhead(lp, i)
            if idx <= m
                k = idx
                get_prim = GLPK.get_row_prim
                get_ub = GLPK.get_row_ub
                get_lb = GLPK.get_row_lb
            else
                k = idx - m
                get_prim = GLPK.get_col_prim
                get_ub = GLPK.get_col_ub
                get_lb = GLPK.get_col_lb
            end

            res = get_prim(lp, k)
            if res > get_ub(lp, k) + eps
                ray[i] = -1
            elseif res < get_lb(lp, k) - eps
                ray[i] = 1
            else
                continue # ray[i] == 0
            end

            if idx <= m
                ray[i] *= GLPK.get_rii(lp, k)
            else
                ray[i] /= GLPK.get_sjj(lp, k)
            end
        end

        GLPK.btran(lp, ray)

        for i = 1:m
            ray[i] /= GLPK.get_rii(lp, i)
        end
    end

    return ray
end

function getunboundedray(lpm::GLPKMathProgModelLP)
    lp = lpm.inner

    if lpm.method == :Simplex || lpm.method == :Exact
    elseif lpm.method == :InteriorPoint
        error("getunboundedray is not available when using the InteriorPoint method")
    else
        error("bug")
    end

    m = GLPK.get_num_rows(lp)
    n = GLPK.get_num_cols(lp)

    ray = zeros(n)

    ur = GLPK.get_unbnd_ray(lp)
    if ur != 0
        if ur <= m
            k = ur
            get_stat = GLPK.get_row_stat
            get_dual = GLPK.get_row_dual
        else
            k = ur - m
            get_stat = GLPK.get_col_stat
            get_dual = GLPK.get_col_dual
            ray[k] = 1
        end

        get_stat(lp, k) != GLPK.BS || error("unbounded ray is dual (use getinfeasibilityray)")

        for (ri, rv) in zip(GLPK.eval_tab_col(lp, ur)...)
            ri > m && (ray[ri - m] = rv)
        end

        if (GLPK.get_obj_dir(lp) == GLPK.MAX) $ (get_dual(lp, k) > 0)
            scale!(ray, -1.0)
        end
    else
        for i = 1:n
            ray[i] = GLPK.get_col_prim(lp, i)
        end
    end

    return ray
end

end


text:
using Documenter
using Jedis

makedocs(
    sitename="Jedis.jl Documentation",
    # format = Documenter.HTML(prettyurls = false),
    pages=[
        "Home" => "index.md",
        "Client" => "client.md",
        "Commands" => "commands.md",
        "Pipelining" => "pipeline.md",
        "Pub/Sub" => "pubsub.md",
        "Locks" => "lock.md"
    ],
    modules=[Jedis]
)

deploydocs(
    repo="github.com/captchanjack/Jedis.jl.git",
    devbranch="main",
    devurl="docs"
)


text:
function segway_rk4(xt::AbstractVector, ut::AbstractVector, dt::Number)
    xt1 = copy(xt)
    segway_rk4!(xt1, ut, dt)
    return xt1
end

function segway_rk4!(xt::AbstractVector, ut::AbstractVector, dt::Number)
    temp = similar(xt)

    #get partial results dy1-dy4
    dy1 = dxdt_segway(xt, ut);
    temp .= dy1.*(dt/2.0) .+ xt;

    dy2 = dxdt_segway(temp, ut);
    temp .= dy2.*(dt/2.0) .+ xt;

    dy3 = dxdt_segway(temp, ut);
    temp .= dy3.*dt .+ xt;

    dy4 = dxdt_segway(temp, ut);

    # update xt
    @. temp = dt*(dy1+2.0*(dy2+dy3)+dy4)/6.0

    #Check if segway has fallen #TODO limits anpassen
    #if -pi/2>=xt[4]+temp[4]
    #    xt[4] = -pi/2
	#	xt[5:7] = 0.0
    #elseif pi/2<=xt[4]+temp[4]
    #    xt[4] = pi/2
	#	xt[5:7] = 0.0
    #else
        # Limit angle to lie between 0 and 2*pi
    @. xt += temp
        #xt[3] = mod(xt[3], 2*pi)
    #end
    return nothing
end

"""
State vector is [x0, y0, phi, alpha, dalpha, v, dphi, ir, il]
"""
function dxdt_segway(x::AbstractVector, u::AbstractVector)

    u1 = u[1];
    u2 = u[2];
    x3 = x[3];
    x4 = x[4];
    x5 = x[5];
    x6 = x[6];
    x7 = x[7];
    t2 = cos(x4);
    t3 = x6.*-2.424242424242424e2;
    t4 = t2.*1.368588515497592e-6;
    t5 = t2.^2;
    t6 = t5.*-1.979293098627009e-7;
    t7 = t4+t6+2.760783558017359e-6;
    t8 = 1.0./t7;
    t9 = x5.*8.0;
    t63 = x7.*1.187878787878788e1;
    t56 = -t63;
    t20 = -t56;
    t10 = -t20;
    t11 = t3+t9+t10;
    t12 = exp(t11);
    t13 = t12+1.0;
    t14 = 1.0./t13;
    t15 = t14.*6.519999999999999e-2;
    t16 = t2.*4.4489247e-4;
    t17 = t16-1.53811157502575e-3;
    t18 = x7.*2.274787878787879e-3;
    t19 = x5.*8.0;
    t21 = t16+3.788361628429753e-4;
    t22 = x5.*6.238865454545455e-3;
    t23 = x6.*-1.890565289256198e-1;
    t24 = cos(x3);
    t25 = sin(x3);
    t26 = sin(x4);
    t27 = x7.^2;
    t28 = x5.^2;
    t61 = x5.*1.532e-3;
    t54 = -t61;
    t31 = -t54;
    t29 = -t31;
    t62 = x6.*-4.642424242424242e-2;
    t55 = -t62;
    t32 = -t55;
    t30 = -t32;
    t33 = x6.*-2.424242424242424e2;
    t34 = u1.*3.3e-2;
    t65 = x7.*9.263769917355372e-3;
    t64 = -t65;
    t39 = -t64;
    t35 = -t39;
    t36 = t22+t23+t34+t35;
    t37 = t16+1.136228822457025e-3;
    t38 = u2.*3.3e-2;
    t40 = t2.*1.368588515497592e-6;
    t41 = t5.*-1.979293098627009e-7;
    t42 = t40+t41+2.760783558017359e-6;
    t43 = 1.0./t42;
    t44 = x5.*-3.064e-3;
    t45 = x6.*9.284848484848484e-2;
    t46 = t26.*1.322543979e-1;
    t53 = x7.*1.187878787878788e1;
    t47 = -t53;
    t48 = t19+t33+t47;
    t49 = exp(t48);
    t50 = t49+1.0;
    t51 = 1.0./t50;
    t52 = t51.*6.519999999999999e-2;
    t57 = t26.^2;
    t58 = t57.*1.1383382179917e-6;
    t59 = t58+4.056771650622809e-6;
    t60 = 1.0./t59;
    t66 = x4.*2.0;
    t67 = sin(t66);
    xDot = [t24.*x6;t25.*x6;x7;x5;t8.*(t15+t44+t45+t46+t27.*t67.*5.2265299265e-4+6.519999999999999e-2./(exp(t3+t19+x7.*1.187878787878788e1)+1.0)-6.519999999999999e-2).*1.916947737868725e-3-t8.*t21.*(t22+t23+t38+x7.*9.263769917355372e-3).*3.819323816679188-t8.*t17.*(t18+x5.*1.532e-3-x6.*4.642424242424242e-2-6.519999999999999e-2./(exp(t3+t19+t20)+1.0)+3.26e-2)+t8.*t17.*(t15+t18+t29+t30-3.26e-2)-t8.*t21.*t36.*3.819323816679188-t8.*t17.*t24.*(t24.*t26.*t27.*4.867e-2+t24.*t26.*t28.*4.867e-2+t25.*x6.*x7+t2.*t25.*x5.*x7.*9.734e-2).*9.141000000000001e-3-t8.*t17.*t25.*(t25.*t26.*t27.*4.867e-2+t25.*t26.*t28.*4.867e-2-t24.*x6.*x7-t2.*t24.*x5.*x7.*9.734e-2).*9.141000000000001e-3;t8.*(t18+t31+t32-6.519999999999999e-2./(exp(t19+t33+x7.*1.187878787878788e1)+1.0)+3.26e-2).*8.825323311693156e-5-t8.*(t18+t29+t30+t52-3.26e-2).*8.825323311693156e-5-t8.*t17.*(t44+t45+t46+t52+6.519999999999999e-2./(exp(t19+t33+t53)+1.0)+t2.*t26.*t27.*1.0453059853e-3-6.519999999999999e-2).*3.3e-2+t8.*t37.*(t22+t23+t38+t39).*1.260376859504132e-1+t8.*t36.*t37.*1.260376859504132e-1+t26.*t27.*t43.*3.926319886687748e-8+t26.*t28.*t43.*3.926319886687748e-8;t60.*(t22+t23+t34+t64).*6.175846611570248e-3-t60.*(t22+t23+t38+t65).*6.175846611570248e-3-t60.*(t18+t54+t55+6.519999999999999e-2./(exp(t3+t9+t56)+1.0)-3.26e-2).*1.617e-3-t60.*(t18+t61+t62-6.519999999999999e-2./(exp(t3+t19+t63)+1.0)+3.26e-2).*1.617e-3-t60.*x7.*(t26.*x6.*1.348159e-2+t67.*x5.*1.0453059853e-3).*1.089e-3];

end

"""
State vector is [x0, y0, phi, alpha, dalpha, v, dphi, ir, il]
"""
function segway_cost(xt::AbstractVector, u::AbstractVector)
    cost = (xt[4]^2 + 0.1*xt[5]^2
            + 10*xt[2]^2 + 10*xt[1]^2 # or squared here...
            + 0.1* xt[7]^2
            + 0.01*(u'u)) # punish for falling over

    reward = 0.0 #reward for moving along the x-axis

    if xt[4]^2 > 2.47 #check if lager than pi/2
        terminal = 1
        cost *= 1.2
    else
        terminal = 0
    end

    return cost-reward, terminal
end

"""
State vector is [x0, y0, phi, alpha, dalpha, v, dphi, ir, il]
"""
function segway_cost2(xt::AbstractVector, u::AbstractVector)
    cost =  10.0*xt[1]^2 + 10.0*xt[2]^2
            + sin(xt[3])^2 + (1.0-cos(xt[3]))^2
            + 10*sin(xt[4])^2 + (1.0-cos(xt[4]))^2
            + 1*xt[5]^2 + 0.1*xt[6]^2
            + 0.1*xt[7]^2 + 0.1*u[1]^2
            + 0.1*u[2]^2

    if xt[4]^2 > 2.47 #check if lager than pi/2
        terminal = 1
        cost *= 10
    else
        terminal = 0
    end

    return cost, terminal
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
# # Molecules
#
# inspired by https://hackage.haskell.org/package/lens-tutorial-1.0.3/docs/Control-Lens-Tutorial.html
#
using Accessors
using Test

molecule = (
    name="water",
    atoms=[
        (name="H", position=(x=0,y=1)), # in reality the angle is about 104deg
        (name="O", position=(x=0,y=0)),
        (name="H", position=(x=1,y=0)),
    ]
)

oc = @optic _.atoms |> Elements() |> _.position.x
res_modify = modify(x->x+1, molecule, oc)

res_macro = @set molecule.atoms |> Elements() |> _.position.x += 1
@test res_macro == res_modify

res_expected = (
    name="water",
    atoms=[
        (name="H", position=(x=1,y=1)),
        (name="O", position=(x=1,y=0)),
        (name="H", position=(x=2,y=0)),
    ]
)

@test res_expected == res_macro

res_set = set(molecule, oc, 4.0)
res_macro = @set molecule.atoms |> Elements() |> _.position.x = 4.0
@test res_macro == res_set

res_expected = (
    name="water",
    atoms=[
        (name="H", position=(x=4.0,y=1)),
        (name="O", position=(x=4.0,y=0)),
        (name="H", position=(x=4.0,y=0)),
    ]
)
@test res_expected == res_set


text:
using Cubature
using Statistics

Bickley(x; U, L, x₀) =  U * sech.((x .- x₀) / L).^2

mutable struct Bickley_Jet{V, P}
 phys::P
    u::V                            # meters/second
    η::V                            # meters
end

"""
Cartesian version.
"""
function Bickley_Jet(geometry::Cartesian;
    phys,
    Uj,
    Ljscale
    )

    Lj = phys.grid.Ly/Ljscale

    LHS(x) = - phys.f₀ / phys.g * Bickley(x; U=Uj, L=Lj, x₀=0)

    u = Bickley(phys.grid.y; U=Uj, L=Lj, x₀=0)
    η = zeros(size(phys.grid.y))
    for i in 1:phys.grid.N
        result = pquadrature(LHS, phys.grid.y[end], phys.grid.y[i])
        η[i] = result[1]
    end
    η = η .- mean(η)

    return Bickley_Jet(phys, u, η)
end

"""
Spherical version.
"""
function Bickley_Jet(geometry::Spherical;
    phys,
    Uj,
    Ljscale
    )

    Lj = phys.grid.Lϕ/Ljscale

    LHS(x) = - (phys.TwoΩ .* sin.(x) .* phys.grid.a .+ Bickley(x; U=Uj, L=Lj, x₀=phys.grid.ϕmid).*tan.(x)) .* Bickley(x; U=Uj, L=Lj, x₀=phys.grid.ϕmid) / phys.g

    u =  Bickley(phys.grid.ϕ; U=Uj, L=Lj, x₀ = phys.grid.ϕmid)
    η = zeros(size(phys.grid.ϕ))
    for i in 1:phys.grid.N
        result = pquadrature(LHS, phys.grid.ϕ[end], phys.grid.ϕ[i])
        η[i] =  result[1]
    end
    η = η .- mean(η)

    return Bickley_Jet(phys, u, η)
end


text:
#import only the needed parts not to have conflict
import NNlib.conv, NNlib.conv!

### NNConv
@doc raw"""
    NNConv(cLayer::CL, Ai::AbstractArray{T,N}) where {T,N, CL <: ConvLayer}

Perform the forward propagation for `cLayer::ConvLayer` using fast implementation of `NNlib`

# Return

- `Dict(:Z => Z, :A => A)`
"""
function NNConv(cLayer::CL, Ai::AbstractArray{T,N}) where {T,N, CL <: ConvLayer}
    padS = paddingSize(cLayer, Ai)
    # axW = axes(cLayer.W)[1:end-2]
    # raxW = reverse.(axW)
    # Z = conv(Ai, cLayer.W[raxW..., :, :], stride = cLayer.s, pad = padS)
    Z = conv(Ai, cLayer.W, stride = cLayer.s, pad = padS, flipped=true)
    axB = axes(cLayer.B)[1:end-2]
    Z .+= cLayer.B[axB..., 1, :]
    actFun = cLayer.actFun
    A = eval(:($actFun($Z)))

    return Dict(:Z => Z, :A => A)
end #function img2colConvolve(cLayer::CL

export NNConv

### dNNConv!

#import only the needed parts not to have conflict
import NNlib.∇conv_data!, NNlib.∇conv_filter, NNlib.DenseConvDims

@doc raw"""
    function dNNConv!(
        cLayer::CL,
        Ai::AbstractArray{T1,N},
        dAi::AbstractArray{T2,N},
        dZ::AbstractArray{T3,N},
    ) where {T1, T2, T3, N, CL <: ConvLayer}

Performs the back propagation for `cLayer::ConvLayer` and save values to the pre-allocated `Array` `dAi` and trainable parameters `W` & `B`

# Arguments

- `cLayer::ConvLayer`

- `Ai::AbstractArray{T1,N}` := the input activation of `cLayer`

- `dAi::AbstractArray{T2,N}` := pre-allocated to hold the derivative of the activation

- `dZ::AbstractArray{T3,N}` := the derivative of the cost to the input of the activation function

# Return

`nothing`

"""
function dNNConv!(
    cLayer::CL,
    Ai::AbstractArray{T1,N},
    dAi::AbstractArray{T2,N},
    dZ::AbstractArray{T3,N},
) where {T1, T2, T3, N, CL <: ConvLayer}

    W = cLayer.W
    padS = paddingSize(cLayer, Ai)
    convdim = DenseConvDims(Ai, W, stride = cLayer.s, padding = padS, flipkernel=true)
    # axW = axes(cLayer.W)[1:end-2]
    # raxW = reverse.(axW)
    # dAi .= ∇conv_data(dZ, W[raxW..., :, :], convdim)
    ∇conv_data!(dAi, dZ, W, convdim)
    cLayer.dW = ∇conv_filter(Ai, dZ, convdim)

    cLayer.dB = sum(permutedims(dZ, [(1:N-2)..., N, N-1]), dims = 1:(N-1))

    return nothing
end #function img2colConvolve(cLayer::Conv1D

export dNNConv!


text:
"""
    exampleobject(nrow, ncol)

Create an example `SummarizedExperiment` object with the specified number of rows and columns.
This is to be used to improve the succinctness of examples and tests.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10)
20x10 SummarizedExperiment
  assays(3): foo bar whee
  rownames: Gene1 Gene2 ... Gene19 Gene20
  rowdata(2): name Type
  colnames: Patient1 Patient2 ... Patient9 Patient10
  coldata(3): name Treatment Response
  metadata(1): version
```
"""
function exampleobject(nrow::Int64, ncol::Int64)
    first_assay = Array{Int64, 2}(undef, (nrow, ncol))
    for i in 1:length(first_assay)
        first_assay[i] = round(rand() * 100)
    end

    second_assay = Array{Float64, 2}(undef, (nrow, ncol))
    for i in 1:length(second_assay)
        second_assay[i] = rand()
    end

    third_assay = Array{Float32, 3}(undef, (nrow, ncol, 2))
    for i in 1:length(third_assay)
        third_assay[i] = -rand()
    end

    assays = OrderedDict{String, AbstractArray}(
        "foo" => first_assay,
        "bar" => second_assay,
        "whee" => third_assay
    )

    tchoices = ["normal", "drug1", "drug2"]
    coldata = DataFrame(
        "name" => [ "Patient" * string(i) for i in 1:ncol ],
        "Treatment" => [ tchoices[Int(floor(rand() * length(tchoices) + 1))] for i in 1:ncol ],
        "Response" => [ rand() for i in 1:ncol ]
    )

    gchoices = ["Protein-coding", "Pseudogene"]
    rowdata = DataFrame(
        "name" => [ "Gene" * string(i) for i in 1:nrow ],
        "Type" => [ gchoices[Int(floor(rand() * length(gchoices) + 1))] for i in 1:nrow ]
    )

    metadata = Dict{String,Any}("version" => "1.1.0")
    return SummarizedExperiment(assays, rowdata, coldata, metadata)
end


text:
@testset "Greedy" begin
    @polyvar x y z w
    f = x^2+y^2+1
    grid = Grid2D((-2, 2), 800)
    A = amoeba(f, grid=grid, alg=Greedy())
    @test any(A.data)

    I = imaginary_projection(f, grid=grid, alg=Greedy())
    @test any(I.data)

    f = x^2+y^2+z^2+1
    F = AmoebaFiber3D(f)
    grid = Grid3D(xlims=(-2, 2), ylims=(-2, 2), zlims=(-2, 2), res=(21, 21, 21))
    A3 = PolynomialAmoebas.amoeba(f, grid=grid, alg=Greedy())
    @test any(A3.data)

    I3 = imaginary_projection(f, grid=grid, alg=Greedy())
    @test any(I3.data)

    f = 1 + z * cis(0.8π) + z^3*w*cis(1.6π) + z*w^2*cis(0.4π) + w*cis(1.2π)
    C = coamoeba(f, alg=Greedy(), resolution=500)
    @test any(C.data)
end


text:
using FastCGI
using Test
using Random
using Sockets

const NLOOPS = 10
const STATS = Dict{String, Any}()

function test_types()
    @testset "utility methods" begin
        @test FastCGI.BYTE(0x0102, 1) === 0x02
        @test FastCGI.BYTE(0x0102, 2) === 0x01
        @test FastCGI.B16(0x02, 0x03) === 0x0203
        @test FastCGI.padding(7) == 1
        @test FastCGI.padding(8) == 0
        @test FastCGI.padding(17) == 7
    end

    @testset "read-write nv lengths" begin
        p = PipeBuffer()
        @test FastCGI._writenvlen(p, UInt8(5)) == 1
        @test FastCGI._readnvlen(p) === UInt32(5)
        @test FastCGI._writenvlen(p, UInt16(127)) == 1
        @test FastCGI._readnvlen(p) === UInt32(127)
        @test FastCGI._writenvlen(p, UInt16(128)) == 4
        @test FastCGI._readnvlen(p) === UInt32(128)
        @test FastCGI._writenvlen(p, UInt16(1024)) == 4
        @test FastCGI._readnvlen(p) === UInt32(1024)
    end

    @testset "read-write types" begin
        p = PipeBuffer()
        reqid = UInt16(1)

        # FCGIBeginRequest
        req = FastCGI.FCGIBeginRequest(FastCGI.FCGIRequestRole.RESPONDER, FastCGI.KEEP_CONN)
        rec = FastCGI.FCGIRecord(FastCGI.FCGIHeaderType.BEGIN_REQUEST, reqid, req)
        @test FastCGI.fcgiwrite(p, rec) == 16
        @test FastCGI.FCGIRecord(p) == rec
        reqid += UInt16(1)

        # FCGIEndRequest
        req = FastCGI.FCGIEndRequest(UInt32(10), FastCGI.FCGIEndRequestProtocolStatus.REQUEST_COMPLETE)
        rec = FastCGI.FCGIRecord(FastCGI.FCGIHeaderType.END_REQUEST, reqid, req)
        @test FastCGI.fcgiwrite(p, rec) == 16
        @test FastCGI.FCGIRecord(p) == rec
        reqid += UInt16(1)

        # FCGIUnknownType
        req = FastCGI.FCGIUnknownType(UInt8(20))
        rec = FastCGI.FCGIRecord(FastCGI.FCGIHeaderType.END_REQUEST, reqid, req)
        @test FastCGI.fcgiwrite(p, rec) == 16
        @test FastCGI.FCGIRecord(p) == rec
        reqid += UInt16(1)

        # FCGIParams
        req = FastCGI.FCGIParams()
        push!(req.nvpairs, FastCGI.FCGINameValuePair("name1", "value1"))
        push!(req.nvpairs, FastCGI.FCGINameValuePair("name2", ""))
        push!(req.nvpairs, FastCGI.FCGINameValuePair("name1", randstring(512)))
        push!(req.nvpairs, FastCGI.FCGINameValuePair(randstring(512), randstring(1024)))
        rec = FastCGI.FCGIRecord(FastCGI.FCGIHeaderType.GET_VALUES, reqid, req)
        @test FastCGI.fcgiwrite(p, rec) == 2096
        @test FastCGI.FCGIRecord(p) == rec
        reqid += UInt16(1)
    end
end

function test_bufferedoutput()
    @testset "buffered output" begin
        out = IOBuffer()
        buffout = FastCGI.BufferedOutput((data)->write(out,data); delay=5.0, bytelimit=10)
        bytes = UInt8[1,2,3,4,5]
        write(buffout, bytes)
        sleep(1)
        @test position(out) == 0
        write(buffout, bytes)
        sleep(1)
        @test position(out) == 10
        write(buffout, bytes)
        sleep(1)
        @test position(out) == 10
        sleep(5)
        @test position(out) == 15
        write(buffout, bytes)
        flush(buffout)
        sleep(1)
        @test position(out) == 20
        close(buffout)
        @test !isopen(buffout.iob)
        sleep(6)
        @test istaskdone(buffout.waiter)
    end
end

function test_clientserver()
    @testset "client-server" begin
        @testset "Unix Domain Socket" begin
            testdir = dirname(@__FILE__)
            socket = joinpath(testdir, "fcgi.socket")

            # start server
            server = FCGIServer(socket)
            @show server
            set_server_runner(FastCGI.ProcessRunner)
            set_server_param("FCGI_MAX_CONNS", "200")
            servertask = @async process(server)
            @test issocket(socket)
            @test isrunning(server)

            # run client
            client = FCGIClient(socket)
            @show client
            headers = Dict("DOCUMENT_ROOT"=>testdir, "SCRIPT_NAME"=>"/hello.sh")
            do_requests(client, headers, "ProcessRunner over Unix Domain Sockets")

            # close
            close(client)
            stop(server)
            @test !isrunning(client)
            @test !isrunning(server)
        end
        @testset "TCP Socket" begin
            testdir = dirname(@__FILE__)
            host = ip"127.0.0.1"
            port = 8989

            # start server
            server = FCGIServer(host, port)
            servertask = @async process(server)
            @test isa(server.lsock, Sockets.TCPServer)
            @test isrunning(server)

            # run client
            client = FCGIClient(host, port)
            headers = Dict("DOCUMENT_ROOT"=>testdir, "SCRIPT_NAME"=>"hello.sh")
            do_requests(client, headers, "ProcessRunner over TCP Sockets")

            # close
            close(client)
            stop(server)
            @test !isrunning(client)
            @test !isrunning(server)
        end
    end
end

function test_clientservertimeout()
    @testset "client-server timeout" begin
        @testset "Unix Domain Socket Timeout" begin
            testdir = dirname(@__FILE__)
            socket = joinpath(testdir, "fcgi.socket")

            # start server
            server = FCGIServer(socket)
            @show server
            set_server_runner(()->FastCGI.ProcessRunner(5))
            servertask = @async process(server)
            @test issocket(socket)
            @test isrunning(server)

            # run client
            client = FCGIClient(socket)
            @show client
            headers = Dict("DOCUMENT_ROOT"=>testdir, "SCRIPT_NAME"=>"/hello.sh")
            do_request(client, headers)
            close(client)
            @test !isrunning(client)

            # run timeout test
            client = FCGIClient(socket)
            @show client
            headers = Dict("DOCUMENT_ROOT"=>testdir, "SCRIPT_NAME"=>"/timeout.sh")
            do_request_timeout(client, headers)
            close(client)
            @test !isrunning(client)

            # close
            stop(server)
            @test !isrunning(server)
        end
        @testset "TCP Socket Timeout" begin
            testdir = dirname(@__FILE__)
            host = ip"127.0.0.1"
            port = 8989

            # start server
            server = FCGIServer(host, port)
            set_server_runner(()->FastCGI.ProcessRunner(5))
            servertask = @async process(server)
            @test isa(server.lsock, Sockets.TCPServer)
            @test isrunning(server)

            # run client
            client = FCGIClient(host, port)
            headers = Dict("DOCUMENT_ROOT"=>testdir, "SCRIPT_NAME"=>"hello.sh")
            do_request(client, headers)
            close(client)
            @test !isrunning(client)

            # run timeout test
            client = FCGIClient(host, port)
            @show client
            headers = Dict("DOCUMENT_ROOT"=>testdir, "SCRIPT_NAME"=>"/timeout.sh")
            do_request_timeout(client, headers)
            close(client)
            @test !isrunning(client)

            # close
            stop(server)
            @test !isrunning(server)
        end
    end
end

function testrunner(params, in, out, err)
    println(out, """Content-type: text/html

<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>Hello World</title>
</head>
<body>
Hello World
<pre>
""")
    for (n,v) in params
        println(out, "$n = $v")
    end
    for (n,v) in ENV
        println(out, "$n = $v")
    end
    println(out, """
</pre>
</body>
</html>
""")

    flush(out)
    close(out)
    close(err)
    return 0
end

function testtimeoutrunner(params, in, out, err)
    sleep(15)
    testrunner(params, in, out, err)
end

function test_functionrunner()
    @testset "Function Runner" begin
        testdir = dirname(@__FILE__)
        socket = joinpath(testdir, "fcgi.socket")

        # start server
        server = FCGIServer(socket)
        @show server
        set_server_runner(FastCGI.FunctionRunner)
        set_server_param("FCGI_MAX_CONNS", "200")
        servertask = @async process(server)
        @test issocket(socket)
        @test isrunning(server)

        # run client
        client = FCGIClient(socket)
        @show client
        headers = Dict("SCRIPT_FILENAME"=>"testrunner")
        do_requests(client, headers, "FunctionRunner over Unix Domain Sockets")

        # close
        close(client)
        stop(server)
        @test !isrunning(client)
        @test !isrunning(server)
    end
end

function do_requests(client, headers, statname)
    STATS[statname] = @timed for idx in 1:NLOOPS
        request = FCGIRequest(; headers=headers, keepconn=true)
        process(client, request)
        @test isrunning(client)
        wait(request.isdone)
        @test isempty(take!(request.err))
        response = take!(request.out)
        @test length(response) > 0
    end
end

function do_request(client, headers)
    request = FCGIRequest(; headers=headers, keepconn=false)
    process(client, request)
    @test isrunning(client)
    wait(request.isdone)
    @test isempty(take!(request.err))
    response = take!(request.out)
    @test length(response) > 0
end

function do_request_timeout(client, headers)
    request = FCGIRequest(; headers=headers, keepconn=false)
    process(client, request)
    @test isrunning(client)
    wait(request.isdone)
    response = take!(request.out)
    @test length(response) == 0
end

function test_functionrunner_timeout()
    @testset "Function Runner Timeout" begin
        testdir = dirname(@__FILE__)
        socket = joinpath(testdir, "fcgi.socket")

        # start server
        server = FCGIServer(socket)
        @show server
        set_server_runner(()->FastCGI.FunctionRunner(5))
        servertask = @async process(server)
        @test issocket(socket)
        @test isrunning(server)

        # run client
        client = FCGIClient(socket)
        @show client
        headers = Dict("SCRIPT_FILENAME"=>"testrunner")
        do_request(client, headers)
        close(client)
        @test !isrunning(client)

        # run timeout test
        client = FCGIClient(socket)
        @show client
        headers = Dict("SCRIPT_FILENAME"=>"testtimeoutrunner")
        do_request_timeout(client, headers)

        close(client)
        @test !isrunning(client)

        # close
        stop(server)
        @test !isrunning(server)
    end
end

test_bufferedoutput()
test_types()
test_clientserver()
test_functionrunner()
test_functionrunner_timeout()
test_clientservertimeout()

println("Times for $NLOOPS requests:")
for (n,v) in STATS
    println(n, ": ", v[2], " sec")
end


text:
function _generate_env_dict(original_env::AbstractDict;
                            JULIA_DEPOT_PATH)
    result = Dict{String, String}()
    if haskey(original_env, "HTTP_PROXY")
        result["HTTP_PROXY"] = deepcopy(original_env["HTTP_PROXY"])
    end
    if haskey(original_env, "HTTPS_PROXY")
        result["HTTPS_PROXY"] = deepcopy(original_env["HTTPS_PROXY"])
    end
    if haskey(original_env, "JULIA_PKG_SERVER")
        result["JULIA_PKG_SERVER"] = deepcopy(original_env["JULIA_PKG_SERVER"])
    end
    result["PATH"] = deepcopy(original_env["PATH"])
    result["JULIA_DEPOT_PATH"] = deepcopy(JULIA_DEPOT_PATH)
    return result
end


text:
module TestFoldsKernelAbstractions
using Test
import CUDA

const TEST_CUDA = CUDA.has_cuda_gpu()
const TEST_GPU = TEST_CUDA

TEST_CUDA && CUDA.allowscalar(false)

find_test(subdir = "") = sort([
    joinpath(subdir, file) for file in readdir(joinpath(@__DIR__, subdir)) if
    match(r"^test_.*\.jl$", file) !== nothing
])

@testset "$file" for file in find_test()
    include(file)
end

@testset "$file" for file in find_test("cpu")
    TEST_GPU && continue
    include(file)
end

@testset "$file" for file in find_test("cuda")
    TEST_CUDA || continue
    include(file)
end

end  # module


text:
module LibNuklear

import Libdl

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Nuklear was not build properly. Please run Pkg.build(\"Nuklear\").")
end
include(depsjl_path)
# Module initialization function
function __init__()
    check_deps()
end

include("CEnum.jl")
using .CEnum

include("ctypes.jl")
export Ctm, Ctime_t, Cclock_t

include(joinpath(@__DIR__, "..", "gen", "libnuklear_common.jl"))
include(joinpath(@__DIR__, "..", "gen", "libnuklear_api.jl"))
include(joinpath(@__DIR__, "..", "gen", "helper_api.jl"))

const nk_false = 0
const nk_true = 1
export nk_false, nk_true

foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "nk_") || startswith(string(s), "NK_")
        @eval export $s
    end
end

end # module


text:
s = """
module Test
export f

function f(x)
    x+1
end

export f #AGAIN!!
end
"""

msgs = lintstr(s)
@test msgs[1].code == :E333
@test contains(msgs[1].message, "duplicate exports of symbol")


text:
# runtests.jl

using Test: @test, @testset, detect_ambiguities
using ImageGeoms

include("helper.jl")

include("core.jl")
include("imresize.jl")
include("makemask.jl")
include("mask.jl")
#include("property.jl")

@testset "ImageGeoms" begin
    @test isempty(detect_ambiguities(ImageGeoms))
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "lttng_tools"
version = v"2.12.4"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://lttng.org/files/lttng-tools/lttng-tools-$(version).tar.bz2",
                  "d729f8c2373a41194f171aeb0da0a9bb35ac181f31afa7e260786d19a500dea1")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/lttng-tools-*
export CPPFLAGS="-I${includedir}"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-man-pages
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("liblttng-ctl", :liblttng_ctl),
    ExecutableProduct("lttng-sessiond", :lttng_sessiond),
    ExecutableProduct("lttng-consumerd", :lttng_consumerd, "lib/lttng/libexec"),
    ExecutableProduct("lttng-crash", :lttng_crash),
    ExecutableProduct("lttng-relayd", :lttng_relayd),
    ExecutableProduct("lttng", :lttng)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="URCU_jll", uuid="aa747835-a391-587f-982f-064ff03f7d29"))
    Dependency(PackageSpec(name="Popt_jll", uuid="e80236cf-ab1d-5f5d-8534-1d1285fe49e8"))
    Dependency(PackageSpec(name="XML2_jll", uuid="02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"))
    Dependency(PackageSpec(name="lttng_ust_jll", uuid="a2826780-45ff-53db-9dda-fd961bc58de1"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Pkg
Pkg.activate(".")
using EcologicalNetworks

N = nz_stream_foodweb()[1]

x = [cascademodel(N) for i in 1:200]


text:
# 这是一个风洞配平程序
# 配平结果要求旋翼没有侧向力
# ——————要求旋翼没有俯仰和滚装力矩
# ——————要求旋翼的桨尖平面与旋翼轴垂直
# 配平结果不要求配平旋翼拉力
# 最后根据一般性的风洞培平要求，决定只通过挥舞量配平操纵量

# 程序首次编辑时间为11/11/2017

function wttrim()
  jacobi = [A0/A_a1 0 A_B1/A_a1;
              γ_B*F0 -B_A1/B_b1  γ_B*F_B1]
  soljaco = inv(jacobi)
  # delthe_lon = 0.1  #纵向周期变距步进
  # delthe_lat = 0.1  #横向周期变距步进

  while abs(Fhub[2])>1e-3 | abs(Mx)>1e-3 | abs(My)>1e-3
    the_lon = the_lon+delthe_lon
    the_lat = the_lat+delthe_lat
    aoaget(v_all,ψ,the_lon,the_lat)
  end

  return the_lon,the_lat

end

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
function mats(p, nx, ny, nu)
	A         = zeros(eltype(p), nx, nx)
	for i = ny:nx-1
		A[i-ny+1,i+1] = 1
	end
	s,e       = 1,nx*ny
	A[:,1:ny] = reshape(p[s:e], nx, ny)
	s,e       = e+1,e+nu*nx
	B         = reshape(p[s:e], nx, nu)
	s,e       = e+1,e+ny*nx
	K         = reshape(p[s:e], nx, ny)
	A,B,K
end
function model_from_params(p, nx, ny, nu)
	A,B,K = mats(p, nx, ny, nu)
	x0    = copy(p[end-nx+1:end])
	sys   = StateSpaceNoise(A,B,K,1.)
	sysf  = SysFilter(sys,x0,similar(x0,ny))
end

function pem_costfun(p,y,u,nx,metric::F) where F # To ensure specialization on metric
	nu,ny = obslength(u),obslength(y)
	model = model_from_params(p, nx, ny, nu)
	return mean(metric(e) for e in prediction_errors(model,y,u))
end
function sem_costfun(p,y,u,nx,metric::F) where F # To ensure specialization on metric
	nu,ny = obslength(u),obslength(y)
	model = model_from_params(p, nx, ny, nu)
	return mean(metric(e) for e in simulation_errors(model,y,u))
end

"""
sys, x0, opt = pem(y, u; nx, kwargs...)

System identification using the prediction error method.

# Arguments:
- `y`: Measurements, either a matrix with time along dim 2, or a vector of vectors
- `u`: Control signals, same structure as `y`
- `nx`: Number of poles in the estimated system. Thus number should be chosen as number of system poles plus number of poles in noise models for measurement noise and load disturbances.
- `focus`: Either `:prediction` or `:simulation`. If `:simulation` is chosen, a two stage problem is solved with prediction focus first, followed by a refinement for simulation focus.
- `metric`: A Function determining how the size of the residuals is measured, default `sse` (e'e), but any Function such as `norm`, `e->sum(abs,e)` or `e -> e'Q*e` could be used.
- `regularizer(p)=0`: function for regularization. The structure of `p` is detailed below
- `solver` Defaults to `Optim.BFGS()`
- `stabilize_predictor=true`: Modifies the estimated Kalman gain `K` in case `A-KC` is not stable by moving all unstable eigenvalues to the unit circle.
- `difficult=false`: If the identification problem appears to be difficult and ends up in a local minimum, set this flag to true to solve an initial global optimization problem to supply a good initial guess. This is expected to take some time.
- `kwargs`: additional keyword arguments are sent to `Optim.Options`.

# Return values
- `sys::StateSpaceNoise`: identified system. Can be converted to `StateSpace` by `convert(StateSpace, sys)`, but this will discard the Kalman gain matrix.
- `x0`: Estimated initial state
- `opt`: Optimization problem structure. Contains info of the result of the optimization problem

## Structure of parameter vecotr `p`
```julia
A = size(nx,ny)
B = size(nx,nu)
K = size(nx,ny)
x0 = size(nx)
p = [A[:];B[:];K[:];x0]
```
"""
function pem(y, u; nx, solver = BFGS(), focus=:prediction, metric=sse, regularizer=p->0, iterations=100, stabilize_predictor=true, difficult=false, kwargs...)
	nu,ny = obslength(u),obslength(y)

	A       = 0.0001randn(nx,ny)
	B       = 0.001randn(nx,nu)
	K       = 0.001randn(nx,ny)
	x0      = 0.001randn(nx)
	p       = [A[:];B[:];K[:];x0]
	options = Options(;iterations=iterations, kwargs...)
	cfp     = p->pem_costfun(p,y,u,nx,metric) + regularizer(p)
	if difficult
		stabilizer = p-> 10000*!stabfun(nx,nu,ny)(p)
		options0 = Options(;iterations=iterations, kwargs...)
		opt = optimize(p->cfp(p)+stabilizer(p), 1000p, ParticleSwarm(n_particles=100length(p)), options0)
		p = minimizer(opt)
	end
	opt     = optimize(cfp, p, solver, options; autodiff = :forward)
	println(opt)
	if focus == :simulation
		@info "Focusing on simulation"
		cf = p->sem_costfun(p,y,u,nx,metric) + regularizer(p)
		opt = optimize(cf, minimizer(opt), NewtonTrustRegion(), Options(;iterations=iterations÷2, kwargs...); autodiff = :forward)
		println(opt)
	end
	model = model_from_params(minimizer(opt), nx, ny, nu)
	if !isstable(model.sys)
		@warn("Estimated system does not have a stable prediction filter (A-KC)")
		if stabilize_predictor
			@info("Stabilizing predictor")
			# model = stabilize(model)
			model = stabilize(model, solver, options, cfp)
		end
	end
	isstable(ss(model.sys)) || @warn("Estimated system is not stable")

	model.sys, copy(model.state), opt
end

function stabilize(model)
	s           = model.sys
	@unpack A,K = s
	C           = s.C
	poles       = eigvals(A-K*C)
	newpoles = map(poles) do p
		ap = abs(p)
		ap <= 1 && (return p)
		p / (ap + sqrt(eps()))
	end
	K2   = ControlSystems.acker(A',C', newpoles)' .|> real
	all(abs(p) <= 1 for p in eigvals(A-K*C)) || @warn("Failed to stabilize predictor")
	s.K .= K2
	model
end

function stabilize(model, solver, options, cfi)
	s = model.sys
	cost = function(p)
		maximum(abs.(eigvals(s.A-p*s.K*s.C)))-0.9999999
	end
	p = fzero(cost,1e-9,1-1e-9)
	s.K .*= p
	model
end

function stabfun(nx,ny,nu)
	function (p)
		model = model_from_params(p,nx,nu,ny)
		isstable(model.sys)
	end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5bf16f26-7191-11eb-1166-73a483cf120a
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)
end

# ╔═╡ 3cfba9ce-7191-11eb-2295-3b4c85656c03
begin
	using MLJ
	using CSV
	using DataFrames
	using PlutoUI
end

# ╔═╡ f3283b68-7191-11eb-3290-33ba4d2aaece
using Random

# ╔═╡ 9dc7a434-7195-11eb-35a1-e9f79dd4ee79
using Test

# ╔═╡ 428fadd2-719a-11eb-3f4b-d98cb77cddbe
using Printf

# ╔═╡ 87aaecb0-7190-11eb-1572-d9110d9ffaaa
md"""
# Week 2: Multiple Regression (Interpretation)
"""

# ╔═╡ c5c92502-7190-11eb-3a75-87b7f156a7ce
md"""
The goal of this first notebook is to explore multiple regression and feature engineering with MLJ functions.

In this notebook we will use data on house sales in King County to predict prices using multiple regression. we will:

  - Use DataFrames to do some feature engineering
  - Use built-in MLJ functions to compute the regression weights (coefficients/parameters)
  - Given the regression weights, predictors and outcome write a function to compute the Residual Sum of Squares (RSS)
  - Look at coefficients and interpret their meanings
  - Evaluate multiple models via RSS

"""

# ╔═╡ 7a50bd00-7191-11eb-2ed6-071e6db68422
md"""
### Load in house sales data
"""

# ╔═╡ 890cb362-7191-11eb-251b-e155ed2d21f2
sales = CSV.File("../../ML_UW_Spec/C02/data/kc_house_test_data.csv";
	header=true) |> DataFrame;

# ╔═╡ b1ae4466-7191-11eb-0c07-29e30cce116a
first(sales, 1)

# ╔═╡ b43046d2-7191-11eb-0978-b1bfe1762756
md"""
### Split data into train/test sets
"""

# ╔═╡ cd45cf32-7191-11eb-1e68-53bccaa9d806
function train_test_split(df; split=0.8, seed=42)
	Random.seed!(seed)
	(nr, nc) = size(df)
	nrp = round(Int, nr * split)
	row_ixes = shuffle(1:nr)
	df_train = view(df[row_ixes, :], 1:nrp, 1:nc)
	df_test = view(df[row_ixes, :], nrp+1:nr, 1:nc)
	(df_train, df_test)
end

# ╔═╡ e054e324-7191-11eb-3e15-c1e0bf73bb0d
sales_train, sales_test = train_test_split(sales);

# ╔═╡ fa6c9bf8-7191-11eb-26b3-1f7f13b91971
size(sales_train), size(sales_test)

# ╔═╡ 0d2ca940-7192-11eb-34ea-ab594b994f41
md"""
### Learning a multiple regression model
"""

# ╔═╡ 14e8024c-7192-11eb-1c76-83a28b6120ee
@load LinearRegressor pkg=MLJLinearModels

# ╔═╡ 72a34b30-7192-11eb-3cda-15a4e9a5244c
describe(sales, :mean, :std, :eltype)

# ╔═╡ c264c3ec-7192-11eb-24ef-c797c42d3671
begin
	example_features = [:sqft_living, :bedrooms, :bathrooms]

	X = select(sales_train, example_features)
	y = sales_train.price

	first(X, 3)
end

# ╔═╡ 1f6d0112-7193-11eb-20a2-efd1550ca8e2
## 3 first house prices
y[1:3]

# ╔═╡ bd39bf8e-7193-11eb-23c0-257408e760a0
size(X), size(y)

# ╔═╡ 45479604-7193-11eb-068e-7ddd263e9be2
## Let's declare a simple multivariate linear regression model:
mdl = MLJLinearModels.LinearRegressor()

# ╔═╡ 89ce6cb2-7193-11eb-25bb-f161f4514a6d
## Let's wrap it in a machine which, in MLJ, is the composition of a model
## and data to apply the model on:

begin
	mach_mr = machine(mdl, X, y)
	fit!(mach_mr)
end

# ╔═╡ d34ab594-7193-11eb-1bed-01a763445a7e
begin
	fp = fitted_params(mach_mr)

	with_terminal() do
		for (name, val) in fp.coefs
    		println("$(rpad(name, 8)):  $(round(val, sigdigits=3))")
		end
		println("Intercept: $(round(fp.intercept, sigdigits=3))")
	end
end

# ╔═╡ 508db812-7194-11eb-3fde-ad654ba107d1
md"""
### Making Predictions
"""

# ╔═╡ 6732d8e0-7194-11eb-3506-93329a94dd84
begin
	ŷ = predict(mach_mr, X)
	(round(ŷ[1], sigdigits=3) , y[1])  # first house predicetd price vs actual price
end

# ╔═╡ c9a0a836-7194-11eb-027a-7bfbfe6cae3c
md"""
### Compute RSS

Now that we can make predictions given the model, let's write a function to compute the RSS [Residual Sum of Squares] of the model.

Let's write a function `get_rss`  to calculate RSS given the machine, data(X), and the outcome (y).
"""

# ╔═╡ 3deaf822-7195-11eb-2bfb-9f51b24982ed
function get_rss(mach, X, y)
    ŷ = predict(mach, X)     # First get the predictions
    diff = y .- ŷ            # Then compute the residuals/errors
    rss = sum(diff .* diff)  # Then square and add them up
    return rss
end

# ╔═╡ a5251374-7195-11eb-368e-5d4814c8a5e2
begin
	X_test = select(sales_test, example_features)
	y_test = sales_test.price

	rss_ex = get_rss(mach_mr, X_test, y_test)
	exp_value = 5.313821155174205e13  # should be 5.313821155174205e13
	@test rss_ex ≈ exp_value
    # println("rss_ex is: $(rss_ex)")
end

# ╔═╡ 83c35d98-7196-11eb-0a75-b952509c4632
md"""
### Create some new features

Let's consider transformations of existing features e.g. the log of the squarefeet or even "interaction" features such as the product of bedrooms and bathrooms.

Create the following 4 new features as column in both test and train data:

  - `bedrooms_squared = bedrooms * bedrooms`
  - `bed_bath_rooms = bedrooms * bathrooms`
  - `log_sqft_living = log(sqft_living)`
  - `lat_plus_long = lat + long`

"""

# ╔═╡ d1cc2b4a-7197-11eb-10fd-8fdadba255d2
begin
	Xtr, Xte = (sales_train, sales_test)
end

# ╔═╡ d0d2312a-7196-11eb-3192-2158f9b7af15
begin
 	nXtr = hcat(Xtr,
		Xtr.bedrooms .* Xtr.bedrooms,
		Xtr.bedrooms .* Xtr.bathrooms,
 		log.(Xtr.sqft_living),
		Xtr.lat .+ Xtr.long;
		makeunique=true);

	rename!(nXtr, :x1 => :bedrooms_squared);
	rename!(nXtr, :x1_1 => :bed_bath_rooms);
	rename!(nXtr, :x1_2 => :log_sqft_living);
	rename!(nXtr, :x1_3 => :lat_plus_long);
end

# ╔═╡ bbe4cfb8-7199-11eb-2c60-edaa2feab075
begin
	nXte = hcat(Xte, Xte.bedrooms .* Xte.bedrooms,
		Xte.bedrooms .* Xte.bathrooms,
		log.(Xte.sqft_living),
		Xte.lat .+ Xte.long;
		makeunique=true);

	rename!(nXte, :x1 => :bedrooms_squared);
	rename!(nXte, :x1_1 => :bed_bath_rooms);
	rename!(nXte, :x1_2 => :log_sqft_living);
	rename!(nXte, :x1_3 => :lat_plus_long);
end

# ╔═╡ fd7ca6fa-7199-11eb-0691-89ea26b44013
md"""

**Quiz Question: What is the mean (arithmetic average) value of your 4 new features on
test data? (round to 2 digits)**

"""

# ╔═╡ 0ac9bda2-719a-11eb-3036-d110c21a2fcc
with_terminal() do
	for attr ∈ [:bed_bath_rooms, :log_sqft_living, :lat_plus_long, :bedrooms_squared]
		@printf("mean of %s is: %3.2f\n", String(attr), mean(nXte[!, attr]))
	end
end

# ╔═╡ afde3d9c-719a-11eb-12cb-8397b9441a65
md"""
### Learning Multiple Models

Now we will learn the weights for three (nested) models for predicting house prices. The first model will have the fewest features the second model will add one more feature and the third will add a few more:
  - Model 1: `squarefeet`, `# bedrooms`, `# bathrooms`, `latitude & longitude`
  - Model 2: add `bedrooms *bathrooms`
  - Model 3: Add `log squarefeet`, `bedrooms squared`, and the (nonsensical) `latitude + longitude`
"""

# ╔═╡ fa0073b6-719a-11eb-3190-dd61c2541ae4
begin
	mod1_features = [:sqft_living, :bedrooms, :bathrooms, :lat, :long];
	mod2_features = [mod1_features..., :bed_bath_rooms];
	mod3_features = [mod2_features..., :bedrooms_squared, :log_sqft_living, :lat_plus_long];
end

# ╔═╡ 5caf09e6-719b-11eb-1652-8938ad992ea3
## model 1
begin
	X1 = select(nXtr, mod1_features)
	mach_mr1 = machine(mdl, X1, y)
	fit!(mach_mr1)

	fp1 = fitted_params(mach_mr1)

	with_terminal() do
		for (name, val) in fp1.coefs
    		println("$(rpad(name, 16)):  $(round(val, sigdigits=3))")
		end
		println("Intercept: $(round(fp1.intercept, sigdigits=3))")
	end
end

# ╔═╡ c56e31be-719b-11eb-022d-f796b95b6b32
## model 2
begin
	X2 = select(nXtr, mod2_features)
	mach_mr2 = machine(mdl, X2, y)
	fit!(mach_mr2)

	fp2 = fitted_params(mach_mr2)

	with_terminal() do
		for (name, val) in fp2.coefs
    		println("$(rpad(name, 16)):  $(round(val, sigdigits=3))")
		end
		println("Intercept: $(round(fp2.intercept, sigdigits=3))")
	end
end

# ╔═╡ d972ae42-719b-11eb-193e-675d543025c1
## model 3
begin
	X3 = select(nXtr, mod3_features)
	mach_mr3 = machine(mdl, X3, y)
	fit!(mach_mr3)

	fp3 = fitted_params(mach_mr3)

	with_terminal() do
		for (name, val) in fp3.coefs
    		println("$(rpad(name, 16)):  $(round(val, sigdigits=3))")
		end
		println("Intercept: $(round(fp3.intercept, sigdigits=3))")
	end
end

# ╔═╡ 266c0a0e-719c-11eb-1692-f74b2a6f5bc4
md"""

**Quiz Question: What is the sign (positive or negative) for the coefficient/weight for 'bathrooms' in model 1?**
  - [x] Positive (+)
  - [ ] Negative (-)

**Quiz Question: What is the sign (positive or negative) for the coefficient/weight for 'bathrooms' in model 2?**
  - [ ] Positive (+)
  - [x] Negative (-)

Think about what this means.

"""

# ╔═╡ 2dfa4aee-719c-11eb-1f89-1babfa6312d2
md"""
### Comparing multiple models

Now that we have learned three models and extracted the model weights we want to evaluate which model is best.

"""

# ╔═╡ 48377010-719c-11eb-19a7-233bf9bb4beb
# function calc_rss(;models=[(:model_1, mach_mr1),
# 			(:model_2, mach_mr2),
# 			(:model_3, mach_mr3)], tdata=sales_train)
# 	for (modn, mach) ∈ models
#     	rss_mod = get_rss(mach, tdata, tdata.price)
#     	@printf("model %10s rss is: %2.7e", String(modn), rss_mod)
# 	end
#   	return
# end

# ╔═╡ 9af9983a-719e-11eb-0a18-f5493f3df6a3
md"""
**Training Data**
"""

# ╔═╡ 06cdc690-719e-11eb-0819-1b5fce1e3a18
begin
	rss_mmr1 = get_rss(mach_mr1, X1, y)
	rss_mmr2 = get_rss(mach_mr2, X2, y)
	rss_mmr3 = get_rss(mach_mr3, X3, y)

	with_terminal() do
		for (ix, rss_) ∈ enumerate((rss_mmr1, rss_mmr2, rss_mmr3))
			@printf("model %2d rss is: %2.7e\n", ix, rss_)
		end
	end
end

# ╔═╡ e60ea370-719f-11eb-3d56-5dcfecf4ab3e
md"""
**Quiz Question: Which model (1, 2 or 3) has lowest RSS on *training* Data?** Is this what you expected?

  - [ ] Model 1
  - [ ] Model 2
  - [x] Model 3

  I would expect this (on training data)
"""

# ╔═╡ 90a5c372-719e-11eb-3483-6140d252b0c0
md"""
**Test Data**
"""

# ╔═╡ b2d5046c-719e-11eb-07db-454247474942
begin
	rss_mmr1_te = get_rss(mach_mr1, select(nXte, mod1_features), sales_test.price)
	rss_mmr2_te = get_rss(mach_mr2, select(nXte, mod2_features), sales_test.price)
	rss_mmr3_te = get_rss(mach_mr3, select(nXte, mod3_features), sales_test.price)

	with_terminal() do
		for (ix, rss_) ∈ enumerate((rss_mmr1_te, rss_mmr2_te, rss_mmr3_te))
			@printf("model %2d rss is: %2.7e\n", ix, rss_)
		end
	end
end

# ╔═╡ fb4038a6-719f-11eb-3e40-cd5aecfb149b
md"""

**Quiz Question: Which model (1, 2 or 3) has lowest RSS on *testing* Data?** Is this what you expected? Think about the features that were added to each model from the previous.

  - [x] Model 1
  - [ ] Model 2
  - [ ] Model 3

"""

# ╔═╡ 0483452a-71a0-11eb-30c0-7d56ed83e837

# ╔═╡ Cell order:
# ╟─87aaecb0-7190-11eb-1572-d9110d9ffaaa
# ╟─c5c92502-7190-11eb-3a75-87b7f156a7ce
# ╠═5bf16f26-7191-11eb-1166-73a483cf120a
# ╠═3cfba9ce-7191-11eb-2295-3b4c85656c03
# ╟─7a50bd00-7191-11eb-2ed6-071e6db68422
# ╠═890cb362-7191-11eb-251b-e155ed2d21f2
# ╠═b1ae4466-7191-11eb-0c07-29e30cce116a
# ╟─b43046d2-7191-11eb-0978-b1bfe1762756
# ╠═f3283b68-7191-11eb-3290-33ba4d2aaece
# ╠═cd45cf32-7191-11eb-1e68-53bccaa9d806
# ╠═e054e324-7191-11eb-3e15-c1e0bf73bb0d
# ╠═fa6c9bf8-7191-11eb-26b3-1f7f13b91971
# ╟─0d2ca940-7192-11eb-34ea-ab594b994f41
# ╠═14e8024c-7192-11eb-1c76-83a28b6120ee
# ╠═72a34b30-7192-11eb-3cda-15a4e9a5244c
# ╠═c264c3ec-7192-11eb-24ef-c797c42d3671
# ╠═1f6d0112-7193-11eb-20a2-efd1550ca8e2
# ╠═bd39bf8e-7193-11eb-23c0-257408e760a0
# ╠═45479604-7193-11eb-068e-7ddd263e9be2
# ╠═89ce6cb2-7193-11eb-25bb-f161f4514a6d
# ╠═d34ab594-7193-11eb-1bed-01a763445a7e
# ╟─508db812-7194-11eb-3fde-ad654ba107d1
# ╠═6732d8e0-7194-11eb-3506-93329a94dd84
# ╟─c9a0a836-7194-11eb-027a-7bfbfe6cae3c
# ╠═3deaf822-7195-11eb-2bfb-9f51b24982ed
# ╠═9dc7a434-7195-11eb-35a1-e9f79dd4ee79
# ╠═a5251374-7195-11eb-368e-5d4814c8a5e2
# ╟─83c35d98-7196-11eb-0a75-b952509c4632
# ╠═d1cc2b4a-7197-11eb-10fd-8fdadba255d2
# ╠═d0d2312a-7196-11eb-3192-2158f9b7af15
# ╠═bbe4cfb8-7199-11eb-2c60-edaa2feab075
# ╟─fd7ca6fa-7199-11eb-0691-89ea26b44013
# ╠═428fadd2-719a-11eb-3f4b-d98cb77cddbe
# ╠═0ac9bda2-719a-11eb-3036-d110c21a2fcc
# ╟─afde3d9c-719a-11eb-12cb-8397b9441a65
# ╠═fa0073b6-719a-11eb-3190-dd61c2541ae4
# ╠═5caf09e6-719b-11eb-1652-8938ad992ea3
# ╠═c56e31be-719b-11eb-022d-f796b95b6b32
# ╠═d972ae42-719b-11eb-193e-675d543025c1
# ╟─266c0a0e-719c-11eb-1692-f74b2a6f5bc4
# ╟─2dfa4aee-719c-11eb-1f89-1babfa6312d2
# ╠═48377010-719c-11eb-19a7-233bf9bb4beb
# ╟─9af9983a-719e-11eb-0a18-f5493f3df6a3
# ╠═06cdc690-719e-11eb-0819-1b5fce1e3a18
# ╟─e60ea370-719f-11eb-3d56-5dcfecf4ab3e
# ╟─90a5c372-719e-11eb-3483-6140d252b0c0
# ╠═b2d5046c-719e-11eb-07db-454247474942
# ╟─fb4038a6-719f-11eb-3e40-cd5aecfb149b
# ╠═0483452a-71a0-11eb-30c0-7d56ed83e837


text:
isabspath(ARGS[1]) || error("Please use an absolute path to the output folder.")
println("Reading experiments from: '$(ARGS[1])'")
data_output_root = ARGS[1]

using OneClassActiveLearning, SVDD
using DataFrames, CSV
using Unmarshal, JSON
using Logging, Statistics
using StatsBase

#### Simplify naming

ALIASES = Dict{Symbol, Dict{String, Union{String, Missing}}}(
    :scale_fct => Dict("min_max_normalize" => "min_max"),
    :combination_fct => Dict("+" => "sum"),
    :init_strategy_scope => Dict("Val{:Subspace}" => "Subspace", "Val{:Global}" => "Global", "" => missing),
    :init_strategy_C => Dict("BoundedTaxErrorEstimate" => "TaxBound",
                             "TaxErrorEstimate" => "Tax",
                             "BinarySearchCStrategy" => "BinarySearch",
                              "FixedCStrategy" => "Fix"),
    :init_strategy_gamma  => Dict("WangGammaStrategy" =>  "Wang",
                                 "RuleOfThumbSilverman" => "Silverman",
                                 "RuleOfThumbScott" => "Scott")
)

function simplify(str, key)
   for (pattern, alias) in ALIASES[key]
       occursin(pattern, str) && return alias
    end
    return str
end

##### Build CSV

array_stats(x) = minimum(x), mean(x), maximum(x)
prefix_named_tuples(x, prefix) = merge(NamedTuple(), [Symbol("$(prefix)_$k")=>v for (k, v) in pairs(x)])

function get_time_mem_stats(r::OneClassActiveLearning.Result)
    if !haskey(r.al_history, :time_fit) || !haskey(r.al_history, :mem_fit) ||
       !haskey(r.al_history, :time_qs) || !haskey(r.al_history, :mem_qs)
        return (time_fit_min=missing, time_fit_mean=missing, time_fit_max=missing,
                mem_fit_min=missing, mem_fit_mean=missing, mem_fit_max=missing,
                time_qs_min=missing, time_qs_mean=missing, time_qs_max=missing,
                mem_qs_min=missing, mem_qs_mean=missing, mem_qs_max=missing)
    end
    time_fit_min, time_fit_mean, time_fit_max = array_stats(values(r.al_history, :time_fit))
    mem_fit_min, mem_fit_mean, mem_fit_max = array_stats(values(r.al_history, :mem_fit))
    time_qs_min, time_qs_mean, time_qs_max = array_stats(values(r.al_history, :time_qs))
    mem_qs_min, mem_qs_mean, mem_qs_max = array_stats(values(r.al_history, :mem_qs))
    return (time_fit_min=time_fit_min, time_fit_mean=time_fit_mean, time_fit_max=time_fit_max,
            mem_fit_min=mem_fit_min, mem_fit_mean=mem_fit_mean, mem_fit_max=mem_fit_max,
            time_qs_min=time_qs_min, time_qs_mean=time_qs_mean, time_qs_max=time_qs_max,
            mem_qs_min=mem_qs_min, mem_qs_mean=mem_qs_mean, mem_qs_max=mem_qs_max)
end

function get_al_summary(r::OneClassActiveLearning.Result, metric::Symbol, metric_short_name::String, n=5, stability_window=10)
    if !haskey(r.al_summary, metric)
        summary = (start_quality=missing, end_quality=missing, maximum=missing, ramp_up=missing, quality_range=missing, total_quality_range=missing,
                   average_end_quality=missing, average_quality_change=missing, average_gain=missing, average_loss=missing,
                   learning_stability=missing, ratio_of_outlier_queries=missing)
    else
        start_quality = r.al_summary[metric][:start_quality]
        end_quality = r.al_summary[metric][:end_quality]
        maximum = r.al_summary[metric][:maximum]
        ramp_up = r.al_summary[metric][:ramp_up][n]
        quality_range = r.al_summary[metric][:quality_range]
        total_quality_range = r.al_summary[metric][:total_quality_range]
        average_end_quality = r.al_summary[metric][:average_end_quality][end-n+1]
        average_quality_change = r.al_summary[metric][:average_quality_change]
        average_gain = r.al_summary[metric][:average_gain]
        average_loss = r.al_summary[metric][:average_loss]
        learning_stability = r.al_summary[metric][:learning_stability][end-stability_window+1]
        ratio_of_outlier_queries = r.al_summary[metric][:ratio_of_outlier_queries]
        summary = (start_quality=start_quality, end_quality=end_quality, maximum=maximum, ramp_up=ramp_up, quality_range=quality_range, total_quality_range=total_quality_range,
                   average_end_quality=average_end_quality, average_quality_change=average_quality_change, average_gain=average_gain, average_loss=average_loss,
                   learning_stability=learning_stability, ratio_of_outlier_queries=ratio_of_outlier_queries)
    end
    return prefix_named_tuples(summary, metric_short_name)
end

function calc_stats(r::OneClassActiveLearning.Result)
    id = string(r.id)
    dir, file_name = splitdir(normpath(r.experiment[:data_file]))
    data_set = r.experiment[:data_set_name]
    split_strategy = r.experiment[:split_strategy_name][1]
    initial_pool_strategy = r.experiment[:initial_pool_strategy_name][1]
    initial_pool_resample_version = r.experiment[:param][:initial_pool_resample_version]
    initial_pool_num_labeled = sum(r.experiment[:param][:initial_pools] .!= "U")
    model = r.experiment[:model][:type]
    if model == "SSAD" && haskey(r.experiment[:model][:param], :κ)
        model = "$(model)_$(r.experiment[:model][:param][:κ])"
    end
    # init_strategy_C = simplify(r.experiment[:model][:init_strategy], :init_strategy_C)
    init_strategy_C = match(r"BoundedTaxErrorEstimate|FixedCStrategy\(\d*.\d*\)", r.experiment[:model][:init_strategy]).match
    init_strategy_gamma = simplify(r.experiment[:model][:init_strategy], :init_strategy_gamma)
    init_strategy_scope = simplify(r.experiment[:model][:init_strategy], :init_strategy_scope)

    qs_raw = split(string(r.experiment[:query_strategy][:type]), ".")[end]
    qs = begin
        if occursin("Sub", qs_raw)
            m = match(r"(\w+)\{(\w+)\}", qs_raw)
            qs_meta = m[:1]
            qs = m[:2]
            qs_scale_fct = simplify(r.experiment[:query_strategy][:param][:scale_fct], :scale_fct)
            qs_combination_fct = simplify(r.experiment[:query_strategy][:param][:combination_fct], :combination_fct)
            (qs=qs, qs_meta=qs_meta, qs_scale_fct=qs_scale_fct, qs_combination_fct=qs_combination_fct)
        else
            (qs=qs_raw, qs_meta=missing, qs_scale_fct=missing, qs_combination_fct=missing)
        end
    end

     subspace_stats = begin
        if !haskey(r.experiment[:model][:param], :subspaces)
            (num_subspaces=missing, max_size_subspaces=missing, min_size_subspaces=missing)
        else
            subspaces = r.experiment[:model][:param][:subspaces]
            (num_subspaces=length(subspaces), max_size_subspaces=maximum(length.(subspaces)), min_size_subspaces=minimum(length.(subspaces)))
        end
    end

    weight_update_strategy = begin
       if haskey(r.experiment[:model][:param], :weight_update_strategy)
           wus = r.experiment[:model][:param][:weight_update_strategy]
           "out_$(wus[:v_Lout])-in_$(wus[:v_Lin])"
       else
           missing
       end
   end

    num_points = r.data_stats.num_observations
    num_dimensions = r.data_stats.num_dimensions
    num_al_iterations = r.experiment[:param][:num_al_iterations]

    exit_code = string(r.status[:exit_code])
    res = (id=id, file_name=file_name, data_set=data_set, subspace_stats..., split_strategy=split_strategy, initial_pool_strategy=initial_pool_strategy,
            initial_pool_resample_version=initial_pool_resample_version, initial_pool_num_labeled=initial_pool_num_labeled, model=model,
            init_strategy_C=init_strategy_C, init_strategy_gamma=init_strategy_gamma, init_strategy_scope=init_strategy_scope, qs...,
            weight_update_strategy=weight_update_strategy,
            num_points=num_points, num_dimensions=num_dimensions, num_al_iterations=num_al_iterations, exit_code=exit_code)
    stats = reduce(merge, [get_time_mem_stats(r),
            get_al_summary(r, :matthews_corr, "m"),
            get_al_summary(r, :cohens_kappa, "c"),
            get_al_summary(r, :auc, "auc"),
            get_al_summary(r, Symbol("auc_fpr_normalized_0.05"), "pauc")])
    return merge(res, stats)
end

function reduce_results(input_path)
    df = Any[]
    for scenario in readdir(input_path)
        if occursin(".csv", scenario)
            continue
        end
        for ds in readdir(joinpath(input_path, scenario, "results"))
            @info ds
            result_path = joinpath(input_path, scenario, "results", ds)
            result_files = joinpath.(result_path, filter(x -> !startswith(x, "."), readdir(result_path)))
            for r in result_files
                @info r
                res = Unmarshal.unmarshal(OneClassActiveLearning.Result, JSON.parsefile(r))
                push!(df,  merge( (scenario=scenario,) , calc_stats(res) ) )
            end
        end
    end
    df_types = mapreduce(x -> [typeof(t) for t in x], (d1, d2) -> [Union{t1,t2} for (t1,t2) in zip(d1, d2)], df)
    df_names = keys(df[1])
    res_df = DataFrame(df_types, [df_names...], 0)
    for row in df
        push!(res_df, row)
    end
    output_file = joinpath(input_path, "$(split(input_path, "/")[end]).csv")
    @info "Writing result to '$(output_file)'"
    CSV.write(output_file, res_df)
    return res_df
end

reduce_results(data_output_root)

max_stars_repo_name: holtri/subsvdd-evaluation

max_issues_repo_name: holtri/subsvdd-evaluation

max_forks_repo_name: holtri/subsvdd-evaluation


text:
cust_z(x,μ,σ) = ismissing(x) ? 0.0 : (x - μ)/σ

function ZScores(x::AbstractArray)
    μ = mean(skipmissing(x))
    σ = std(skipmissing(x))
    cust_z.(x,μ,σ)
end

delta(vec) = vec - lag(vec)

function repeat_by!(df,f,b,cc, n)
    c_list = [Symbol(join([n, x],"_")) for x in cc]
    for c in c_list
        df[!,c] = Vector{Union{Float64,Missing}}(undef,nrow(df))
    end
    combine(groupby(df,b)) do dd
        for (col,rcol) in zip(cc,c_list)
            dd[:,rcol] = f(dd[:,col])
        end
    end
    return df
end

function elaborate(df)
    # df = get_stats()
    repeat_by!(df,delta,:name,deltas_cols,"delta")
    cc = deltas_cols
    c_list = [Symbol(join(["delta", x],"_")) for x in deltas_cols]
    cc = vcat(c_list,z_cols)
    repeat_by!(df,ZScores,:date,cc,"zscore")
    df.zscore_delta_arena_place = df.zscore_delta_arena_place .* -1
    df.zscore_total = df.zscore_delta_arena_victories .+
        # df.zscore_delta_arena_place .+
        # df.zscore_delta_arena_defences .+
        df.zscore_delta_maxed_heroes .+
        df.zscore_clan_tokens .+
        df.zscore_raid_points .+
        df.zscore_portal_stones
    cc = names(df)
    filter!(t -> occursin("zscore",string(t)),cc)
    repeat_by!(df,t->tiedrank(t;rev = true),:date,cc,"rank")
    return df
end

function update_results(clan::AbstractDict)
    full = get_stats(clan)
    elab = elaborate(full)
    return elab
end


text:
# Setup the paths -
const _PATH_TO_ROOT = pwd()
const _PATH_TO_SRC = "$(_PATH_TO_ROOT)/src"
const _PATH_TO_CONFIG = "$(_PATH_TO_ROOT)/config"

# endpoint -
const _TEXT_MAGIC_REST_ENDPOINT = "https://rest.textmagic.com/api/v2"

# load external modules -
using HTTP
using DataFrames
using CSV
using TOML
using JSON
using FileWatching
using Logging

# load my codes -
include("./Types.jl")
include("./Files.jl")
include("./Factory.jl")
include("./Network.jl")


text: const CreateOption = String


text:
function get_chr_sizes_dict(chrom_sizes_path)
    chrom_sizes=readdlm(chrom_sizes_path)
    chrom_sizes_dict = Dict{ASCIIString,Int64}()
    size(chrom_sizes,1)
    for i=1:size(chrom_sizes,1)
        chrom_sizes_dict[ chrom_sizes[i,1] ] = chrom_sizes[i,2]
    end
    return chrom_sizes_dict
end


text:
print_version() = println(COSMO.version())

function print_header(ws::COSMO.Workspace{T}) where {T <: AbstractFloat}

	m, n = ws.p.model_size
	settings = ws.settings
	settings.scaling > 0 ? scaling_status = "on" : scaling_status = "off"

	println("-"^66 * "\n" * " "^10 * "COSMO v0.8.1 - A Quadratic Objective Conic Solver\n" * " "^25 * "Michael Garstka\n"  * " "^16 * "University of Oxford, 2017 - 2021\n" * "-"^66 * "\n")

	println("Problem:  x ∈ R^{$(n)},\n          constraints: A ∈ R^{$(m)x$(n)} ($(count(!iszero, ws.p.A)) nnz),\n          matrix size to factor: $(n + m)x$(n + m),\n          Floating-point precision: $(T)")
	for (iii, set) in enumerate(sort(ws.p.C.sets, by = x -> -x.dim))
		if iii == 1
			print("Sets:"*" "^5)
			print_set(set)
		else
			print(" "^10)
			print_set(set)
		end
		if iii == 5
			println(" "^10*"... and $(length(ws.p.C.sets)-5) more")
			break
		end
	end
	if ws.ci.decompose && ws.ci.num_decomposable > 0
		println("Decomp:   Num of original PSD cones: $(ws.ci.num_psd_cones)\n" * " "^10*"Num decomposable PSD cones: $(ws.ci.num_decomposable)\n" * " "^10*"Num PSD cones after decomposition: $(ws.ci.num_decom_psd_cones)\n" * " "^10*"Merge Strategy: $(stringify(ws.settings.merge_strategy))")
	end
	println("Settings: ϵ_abs = $(@sprintf("%.1e", settings.eps_abs)), ϵ_rel = $(@sprintf("%.1e", settings.eps_rel)),\n" * " "^10 * "ϵ_prim_inf = $(@sprintf("%.1e", settings.eps_prim_inf)), ϵ_dual_inf = $(@sprintf("%.1e", settings.eps_dual_inf)),\n" * " "^10 * "ρ = $(@sprintf("%.3g", settings.rho)), σ = $(@sprintf("%.3g", settings.sigma)), α = $(@sprintf("%.3g", settings.alpha)),\n" * " "^10 * "max_iter = $(settings.max_iter),\n" * " "^10 * "scaling iter = $(settings.scaling) ($(scaling_status)),\n" * " "^10 * "check termination every $(settings.check_termination) iter,\n" * " "^10 * "check infeasibility every $(settings.check_infeasibility) iter,\n" *	 " "^10 * "KKT system solver: $(print_lin_sys(settings.kkt_solver.ObjectType))")

	print_accelerator(ws.accelerator, ws.settings.safeguard, ws.settings.safeguard_tol, tab = 10)

	println("Setup Time: $(round.(ws.times.setup_time*1000; digits=2))ms\n")
	println("Iter:\tObjective:\tPrimal Res:\tDual Res:\tRho:")
	nothing
end

function print_iteration(ws::COSMO.Workspace, iter::Int, cost::T, r_prim::T, r_dual::T) where {T <: AbstractFloat}
	settings = ws.settings
	if mod(iter, 1) == 0 || iter == 1 || iter == 2 || iter == settings.max_iter
		if mod(iter, settings.check_termination) == 0 || iter == 1
			@printf("%d\t%.4e\t%.4e\t%.4e\t%.4e\n", iter, cost, r_prim, r_dual, ws.ρ)
		else
			@printf("%d\t%.4e\t ---\t\t\t---\n", iter, cost)
		end
	end
	nothing
end

function stringify(merge_strategy::Union{Type{<: AbstractMergeStrategy}, OptionsFactory{<: AbstractMergeStrategy}})
	if typeof(merge_strategy) <: OptionsFactory
		return stringify(merge_strategy.ObjectType)
	else
		s = string(merge_strategy)
		# remove the COSMO prefix
		return s[7:end]
	end
end

function print_result(status::Symbol, iter::Int, safeguarding_iter::Int, cost::T, rt::Float64, safeguard::Bool) where {T <: AbstractFloat}
	print("\n" * "-"^66 * "\n>>> Results\nStatus: ")
	if status == :Solved
		result_color = :green
	else
		result_color = :red
	end
	printstyled("$(status)\n", color = result_color)
	print("Iterations: $(iter)")
	if safeguard && safeguarding_iter > 0
		println(" (incl. $(safeguarding_iter) safeguarding iter)")
	else
		print("\n")
	end

	println("Optimal objective: $(@sprintf("%.4g", cost))\nRuntime: $(round.(rt; digits = 3))s ($(round.(rt * 1000; digits = 2))ms)\n")
	nothing
end

function print_lin_sys(s::Type{<:AbstractKKTSolver})
	lin_sys_dict = Dict("QdldlKKTSolver" => "QDLDL", "CholmodKKTSolver" => "CHOLMOD", "PardisoDirectKKTSolver" => "Pardiso (direct)", "PardisoIndirectKKTSolver" => "Pardiso (indirect)",  "MKLPardisoKKTSolver" => "MKL Pardiso")
	return get(lin_sys_dict, string(s), string(s))
end

"Print information about the accelerator at the start."
function print_accelerator(s::CA.AndersonAccelerator{T, BT, ME, RE}, safeguard::Bool, safeguard_tol; tab::Int64) where {T, BT, ME, RE}
		me = ME == CA.RestartedMemory ? "RestartedMemory" : "RollingMemory"
		if BT == CA.Type1
			bt = "Type1"
		elseif BT == CA.Type2{NormalEquations}
			bt = "Type2{NormalEquations}"
		elseif BT == CA.Type2{QRDecomp}
			bt = "Type2{QRDecomp}"
		end

		println("Acc:" * " "^(tab - 4)  * "Anderson $(bt),\n" * " "^tab * "Memory size = $(s.mem), $(me),	\n" * " "^tab * "Safeguarded: $(safeguard), tol: $(safeguard_tol)")
end

print_accelerator(s::CA.AbstractAccelerator, safeguarded::Bool, safeguarding_tol; tab::Int64) = println("Acc:" * " "^(tab - 4) * "Unknown Accelerator")
print_accelerator(s::CA.EmptyAccelerator,  safeguarded::Bool, safeguarding_tol; tab::Int64) = println("Acc:" * " "^(tab - 4) * "no acceleration")

function print_set(set::COSMO.AbstractConvexSet)
	set_name = split(string(typeof(set).name), ".")[end]
	println("$(set_name) of dim: $(set.dim)")
end

function print_set(set::Union{PsdCone{T}, DensePsdCone{T}, PsdConeTriangle{T}, DensePsdConeTriangle{T}}) where {T <: AbstractFloat}
	set_name = split(string(typeof(set).name), ".")[end]
	N = set.sqrt_dim
	println("$(set_name) of dim: $(set.dim) ($(N)x$(N))")
end