text:
## -------------------------------------------------------------------
# Xexp Xmax
function plot_Xexp_Xmax_study()
    methodL = :FBA_Z_FIX_MIN_VG_MIN_COST
    methodU = :FBA_Z_FIX_MAX_VG_MIN_COST

    p = plot(;
        xlabel = _textbf("exp. X (log scale)"),
        ylabel = _textbf("theo. max. X (log scale)"),
    )

    all_vals = []
    for (src_name, iJR, Data) in [
        ("Kayser", ChK.iJR904, ChK.KayserData),
        ("Nanchen", ChN.iJR904, ChN.NanchenData),
        ("Folsom", ChF.iJR904, ChF.FolsomData),
    ]
        src = nameof(Data)
        marker = (9, source_markers[Data])

        @info("Doing", src)

        DAT = ChE.load_DAT(src)
        EXPS = DAT[:EXPS]
        ug_mins = DAT[methodL, :flx, "GLC", EXPS]

        Xexps = [Data.val(:X, exp) for exp in EXPS]
        cgs = [Data.cval(:GLC, exp) for exp in EXPS]
        Ds = [Data.val(:D, exp) for exp in EXPS]
        Xmaxs = ((cgs .* Ds) ./ ug_mins) .|> abs

        Xexps = log10.(Xexps)
        Xmaxs = log10.(Xmaxs)

        scatter!(p, Xexps, Xmaxs;
            label = _textbf(src_name),
            color = :black,
            marker, alpha = 0.8,
        )
        # # fake label for spacing
        # scatter!(p, [Xexps[end]], [Xmaxs[end]]; label = " ", alpha = 0.0)

        plot!(p; legend = :topleft)

        push!(all_vals, Xexps...)
    end

    sort!(all_vals)
    plot!(p, all_vals, all_vals;
        label = "", lw = 3, ls = :dash,
        alpha = 0.7, color = :black
    )
    plot!(;legend = :bottomright)
    # sfig(ChE, p,
    #     "Xexp_Xmax_study", ".png"
    # )

    return p
end


text:
"""
    compute_frames!(actions::Vector{AA}; last_frames=nothing) where AA <: AbstractAction

Set action.frames.frames to the computed frames for each action in actions.
"""
function compute_frames!(
    actions::Vector{AA};
    last_frames = nothing,
) where {AA<:AbstractAction}
    for action in actions
        if last_frames === nothing && get_frames(action) === nothing
            throw(ArgumentError("Frames need to be defined explicitly in the initial
            AbstractAction like Action/BackgroundAction or SubAction."))
        end
        if get_frames(action) === nothing
            set_frames!(action, last_frames)
        end
        last_frames = get_frames(action)
    end
end

"""
    get_current_setting()

Return the current setting of the current action
"""
function get_current_setting()
    action = CURRENT_ACTION[1]
    return action.current_setting
end

"""
    get_interpolation(frames::UnitRange, frame)

Return a value between 0 and 1 which represents the relative `frame` inside `frames`.
"""
function get_interpolation(frames::UnitRange, frame)
    frame == last(frames) && return 1.0
    t = (frame - first(frames)) / (length(frames) - 1)
    # makes sense to only allow 0 ≤ t ≤ 1
    t = min(1.0, t)
end

"""
    get_interpolation(action::AbstractAction, frame)

Return the value of the `action.anim` Animation based on the relative frame given by
`get_interpolation(get_frames(action), frame)`
"""
function get_interpolation(action::AbstractAction, frame)
    t = get_interpolation(get_frames(action), frame)
    if !(action.anim.frames[end].t ≈ 1)
        @warn "Animations should be defined from 0.0 to 1.0"
    end
    return at(action.anim, t)
end

function isapprox_discrete(val; atol = 1e-4)
    return isapprox(val, round(val); atol = atol)
end

ext: jl


text:
using Test

import DivPopGa.ClusteredGa as CLGA
import DivPopGa.Summary as Summary

@testset "Summary" begin
    function costfn(vals)
        return abs(vals[1] - 3.141592) + abs(vals[2] - 2.71828)
    end

    crossfn = CLGA.makelinearcrossover(costfn)
    mutatefn = CLGA.makenormalmutation(10.0,  #stddev
                                      0.10   #mutation prob
                )
    result = CLGA.ga(
        100, # popsize
        500, # generations
        [-100.0, -100.0], #lower
        [100.0, 100.0],   #upper
         costfn, #costfunction
         crossfn, #crossoverfunction
         mutatefn, #mutation function
         CLGA.GA_TYPE_CLASSIC # classical selection
    )

    s = Summary.summary(result)

    @test s isa Dict
    @test s["min_cost"] <= s["q25_cost"]
    @test s["q25_cost"] <= s["median_cost"]
    @test s["median_cost"] <= s["q75_cost"]
    @test s["q75_cost"] <= s["max_cost"]
    @test s["min_cost"] <= s["mean_cost"] <= s["max_cost"]
    @test s["best_solution"] isa CLGA.Chromosome
end

max_stars_repo_name: akadal/divpopga

max_issues_repo_name: akadal/divpopga

max_forks_repo_name: akadal/divpopga


text:
using LinearAlgebra
using Random
using StatsBase
using Combinatorics
import Base.print_matrix

function force(M, i)
    f = zeros(size(M, 1))
    for j in 1:K
        if i == j
            continue
        end
        #d = M[:, i] - M[:, j] .+ 0.00000000000001#+ (rand(N)* 2.0 .- 1.0)*0.00000001 # + rand(N)*0.025
        #d = M[:, i] - M[:, j] + rand(N)*0.025
        d = M[:, i] - M[:, j] + (rand(N) * 2.0 .- 1.0)*0.0125
        f += d * 1/(norm(d)^4.5)
    end
    return f
end

function run(N, C, R, K, ITER_NUM)
    # initial set of tickets(randomly assigned)
    M = zeros(N, K)
    for i in 1:K
        M[randperm(N)[1:C], i] .= 1
    end

    center = ones(N) * C/N

    #coverage = sum(M, dims=2)
    combination_indices = collect(combinations(1:C, R))
    coverage = zeros([N for i in 1:R]...)
    for col in eachcol(M)
        numbers = sort!(findall(!iszero, col))
        for i in combination_indices
            coverage[numbers[i]...] += 1.0
        end
    end

    #return [findall(!iszero, r) for r in eachcol(M)]
    for n in 1:ITER_NUM
        F = zeros(N, K)
        for i in 1:K
            F[:, i] = force(M, i)
        end
        f_mags = [norm(col) for col in eachcol(F)]
        #i = sample(Weights(f_mag))
        i = argmax(f_mags)
        d = normalize(M[:, i] + F[:, i] - center) + center - M[:, i]
        #d_max_index = sample(Weights(d))
        #d_min_index = sample(Weights(-d))
        d_max_index = argmax(d)
        d_min_index = argmin(d)
        if d_max_index != d_min_index && M[d_max_index, i] != 1.0
            M[d_max_index, i] = 1.0
            M[d_min_index, i] = 0.0
            new_f_mag = norm(force(M, i))
            if new_f_mag >= f_mags[i]
                #println(n, "th undo force: ", f_mags[i], "->", new_f_mag)
                M[d_max_index, i] = 0.0
                M[d_min_index, i] = 1.0
            end
        end
        #=
        if d_max_index != d_min_index && M[d_max_index, i] != 1.0
            last_nh = count(x -> x >= 1, coverage)
            last_numbers = sort!(findall(!iszero, vec(M[:, i])))
            for k in combination_indices
                coverage[last_numbers[k]...] -= 1.0
            end
            M[d_max_index, i] = 1.0
            M[d_min_index, i] = 0.0
            new_numbers = sort!(findall(!iszero, vec(M[:, i])))
            for k in combination_indices
                coverage[new_numbers[k]...] += 1.0
            end
            new_nh = count(x -> x >= 1, coverage)
            #println(n, "th hueristic: ", last_nh, "->", new_nh)
            if new_nh < last_nh
                #println("undo..", M[d_max_index, i], M[d_min_index, i])
                M[d_max_index, i] = 0.0
                M[d_min_index, i] = 1.0
                for k in combination_indices
                    coverage[new_numbers[k]...] -= 1.0
                    coverage[last_numbers[k]...] += 1.0
                end
            end
        end=#
    end
    return [findall(!iszero, r) for r in eachcol(M)]
end

N = parse(Int64, ARGS[1])
C = parse(Int64, ARGS[2])
K = parse(Int64, ARGS[3])
R = parse(Int64, ARGS[4])
ITER_NUM = parse(Int64, ARGS[5])
SEED = parse(Int64, ARGS[6])

#const N = 45
#const C = 6
#const R = 3
#const K = 50
#const SEED = 1
#const ITER_NUM = 1000
Random.seed!(SEED)

res = run(N, C, R, K, ITER_NUM)
for s in res
    for num in s
        print(num, ' ')
    end
    println()
end


text:
function _testset_output_dim(scale)
    n = 20
    data = collect(range(-1, 1; length=n))
    û = tvdiff(data, 1, 0.2; scale=scale)

    @testset "Dimensions" begin
        @test length(û) == n # output dim should equal input dim
    end
end


text:
# Autogenerated wrapper script for GSL_jll for i686-w64-mingw32
export gsl_histogram, gsl_randist, libgsl, libgslcblas

JLLWrappers.@generate_wrapper_header("GSL")
JLLWrappers.@declare_executable_product(gsl_histogram)
JLLWrappers.@declare_executable_product(gsl_randist)
JLLWrappers.@declare_library_product(libgsl, "libgsl-25.dll")
JLLWrappers.@declare_library_product(libgslcblas, "libgslcblas-0.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        gsl_histogram,
        "bin\\gsl-histogram.exe",
    )

    JLLWrappers.@init_executable_product(
        gsl_randist,
        "bin\\gsl-randist.exe",
    )

    JLLWrappers.@init_library_product(
        libgsl,
        "bin\\libgsl-25.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgslcblas,
        "bin\\libgslcblas-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file is a part of UpROOT.jl, licensed under the MIT License (MIT).

"""
    TTree <: AbstractVector{Any}

`UpROOT.TTree` is a wrapper around Python objects with mix-in
`uproot.tree.TTreeMethods`. It behaves like a Julia `AbstractVector` and
`Tables.Table` (with column access).

Limitations: Write access is not implemented yet.
"""
struct TTree <: AbstractVector{NamedTuple}
    pyobj::PyObject

    function TTree(x::PyObject)
        if !pybuiltin(:isinstance)(x, uproot.tree.TTreeMethods)
            throw(ArgumentError("Python object $x is not a ROOT TTree"))
        end
        new(x)
    end
end

export TTree

pyobj(tree::TTree) = Base.getfield(tree, :pyobj)

Base.propertynames(tree::TTree) = Symbol.(pyobj(tree).keys())

Base.getproperty(tree::TTree, name::Symbol) = TBranch(get(pyobj(tree), String(name)))

Base.size(tree::TTree) = (pyobj(tree).numentries,)

Base.IndexStyle(tree::TTree) = IndexLinear()

function Base.getindex(tree::TTree, idxs::AbstractUnitRange)
    @boundscheck checkbounds(tree, idxs)
    cols = pyobj(tree).arrays(entrystart = first(idxs) - 1, entrystop = last(idxs))
    Table(_dict2nt(cols))
end

Base.getindex(tree::TTree, i::Integer) = first(tree[i:i])

Base.getindex(tree::TTree, ::Colon) = getindex(tree, eachindex(tree))

Tables.istable(::Type{TTree}) = true

Tables.columnaccess(::Type{TTree}) = true

Tables.columns(tree::TTree) = tree # Dict(broadcast(item -> Symbol(item[1]) => TBranch(item[2]), pyobj(tree).items()))

function Tables.schema(tree::TTree)
    props = Base.propertynames(tree)
    Tables.Schema(props, map(_ -> Any, props))
end

"""
    TBranch <: AbstractVector{Any}

`UpROOT.TBranch` is a wrapper around Python objects with mix-in
`uproot.tree.TBranchMethods`. It behaves like a Julia `AbstractVector`. Ff
the branch has children, it also behaves like a `Tables.Table` (with column
access).

Limitations: Write access is not implemented yet.
"""
struct TBranch <: AbstractVector{Any}
    pyobj::PyObject

    function TBranch(x::PyObject)
        if !pybuiltin(:isinstance)(x, uproot.tree.TBranchMethods)
            throw(ArgumentError("Python object $x is not a ROOT TBranch"))
        end
        new(x)
    end
end

export TBranch

Base.propertynames(branch::TBranch) = Symbol.(pyobj(branch).keys())

Base.getproperty(branch::TBranch, name::Symbol) = TBranch(get(pyobj(branch), String(name)))

Base.size(branch::TBranch) = (pyobj(branch).numentries,)

Base.IndexStyle(branch::TBranch) = IndexLinear()

function Base.getindex(branch::TBranch, idxs::AbstractUnitRange)
    @boundscheck checkbounds(branch, idxs)
    py2jl(pyobj(branch).array(entrystart = first(idxs) - 1, entrystop = last(idxs)))
end

Base.getindex(branch::TBranch, i::Integer) = first(branch[i:i])

Base.getindex(branch::TBranch, ::Colon) = getindex(branch, eachindex(branch))

Tables.istable(::Type{TBranch}) = true

Tables.columnaccess(::Type{TBranch}) = true

Tables.columns(branch::TBranch) = branch # Dict(broadcast(item -> Symbol(item[1]) => TBranch(item[2]), pyobj(branch).items()))

function Tables.schema(branch::TBranch)
    props = Base.propertynames(branch)
    Tables.Schema(props, map(_ -> Any, props))
end

function read_ttree(tree::PyObject, branchnames::Array{<:AbstractString})
    tree_data = tree.arrays(branchnames)

    df = DataFrame()
    for bn in branchnames
        df[Symbol(bn)] = py2jl(tree_data[bn])
    end

    df
end

function read_ttree(filename::AbstractString, treename::AbstractString)
    file = uproot.open(filename)
    tree = get(file, treename)
    branchnames = tree.keys()
    read_ttree(tree, branchnames)
end

function read_ttree(filename::AbstractString, treename::AbstractString, branchnames::Array{<:AbstractString})
    file = uproot.open(filename)
    tree = get(file, treename)
    read_ttree(tree, branchnames)
end

function read_ttree(filename::AbstractString, treename::AbstractString, branchnames::Array{Symbol})
    read_ttree(filename, treename, String.(branchnames))
end


text:
module Powersense

using LinearAlgebra
using SparseArrays
using Printf
using Logging
import MathOptInterface, JuMP, PowerModels

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

export create_PowersenseData, display_preprocess_info, create_opf_model!, run_opf!

include("opt/status.jl")
include("opt/parameters.jl")
include("opt/model.jl")
include("opt/algorithms.jl")
include("opt/MOI_wrapper.jl")

include("opf/PowersenseData.jl");
include("opf/formulations.jl");
end


text:
# Using the Laplace Type Estimator formualation of MSM
# See: https://github.com/JulienPascal/MCMC_Approach_Classical_Estimation/tree/master
#-------------------------------------------------------------------------------
using ClusterManagers
using Distributed

OnCluster = false #set to false to run locally
addWorkers = true #set to false to run serially
println("OnCluster = $(OnCluster)")

# Current number of workers
#--------------------------
currentWorkers = nworkers()
println("Initial number of workers = $(currentWorkers)")

# Increase the number of workers available
#-----------------------------------------
maxNumberWorkers = 10
if addWorkers == true
	if OnCluster == true
	  addprocs(SlurmManager(maxNumberWorkers))
	else
	  addprocs(maxNumberWorkers)
	end
end

# Sanity checks
#-------------
hosts = []
pids = []
for i in workers()
	host, pid = fetch(@spawnat i (gethostname(), getpid()))
	println("Hello I am worker $(i), my host is $(host)")
	push!(hosts, host)
	push!(pids, pid)
end

currentWorkers = nworkers()
println("Number of workers = $(currentWorkers)")

using Plots
using ParallelDataTransfer
using DataFrames
using RobustPmap #Necessary for AffineInvariantMCMC
using CSV
@everywhere using MSM
@everywhere using DataStructures
@everywhere using OrderedCollections
@everywhere using Distributions
@everywhere using Random
@everywhere using DataStructures
@everywhere using Statistics
@everywhere using LinearAlgebra
@everywhere using AffineInvariantMCMC
@everywhere using LinearAlgebra

# Generate simulated data
Random.seed!(1234)  #for replicability reasons
T = 100000          #number of periods
P = 2               #number of dependent variables
beta0 = rand(P)     #choose true coefficients by drawing from a uniform distribution on [0,1]
alpha0 = rand(1)[]  #intercept
theta0 = 0.0        #coefficient to create serial correlation in the error terms
println("True intercept = $(alpha0)")
println("True coefficient beta0 = $(beta0)")
println("Serial correlation coefficient theta0 = $(theta0)")
true_vals = [alpha0, beta0[1], beta0[2]]
# Generation of error terms
# row = individual dimension
# column = time dimension
U = zeros(T)
d = Normal()
U[1] = rand(d, 1)[] #first error term
# loop over time periods
for t = 2:T
    U[t] = rand(d, 1)[] + theta0*U[t-1]
end

# Let's simulate the dependent variables x_t
x = zeros(T, P)

d = Uniform(0, 5)
for p = 1:P
    x[:,p] = rand(d, T)
end

# Let's calculate the resulting y_t
y = zeros(T)

for t=1:T
    y[t] = alpha0 + x[t,1]*beta0[1] + x[t,2]*beta0[2] + U[t]
end

# Send simulated data to workers
sendto(workers(), y=y)
sendto(workers(), x=x)

# Visualize data
p1 = scatter(x[1:100,1], y[1:100], xlabel = "x1", ylabel = "y", legend=:none, smooth=true)
p2 = scatter(x[1:100,2], y[1:100], xlabel = "x2", ylabel = "y", legend=:none, smooth=true)
plot(p1, p2)

# Define locally
optionsSMM = MSMOptions(maxFuncEvals=1000, globalOptimizer = :dxnes, localOptimizer = :NelderMead)
myProblem = MSMProblem(options = optionsSMM);

# Send to workers
sendto(workers(), optionsSMM=optionsSMM)
sendto(workers(), myProblem=myProblem)

# Priors
dictPriors = OrderedDict{String,Array{Float64,1}}()
dictPriors["alpha"] = [0.5, 0.001, 1.0]
dictPriors["beta1"] = [0.5, 0.001, 1.0]
dictPriors["beta2"] = [0.5, 0.001, 1.0]

# Empirical moments
dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
dictEmpiricalMoments["mean"] = [mean(y)] #informative on the intercept
dictEmpiricalMoments["mean^2"] = [mean(y.^2)] #informative on the intercept
dictEmpiricalMoments["mean_x1y"] = [mean(x[:,1] .* y)] #informative on betas
dictEmpiricalMoments["mean_x2y"] = [mean(x[:,2] .* y)] #informative on betas
dictEmpiricalMoments["mean_x1y^2"] = [mean((x[:,1] .* y).^2)] #informative on betas
dictEmpiricalMoments["mean_x2y^2"] = [mean((x[:,2] .* y).^2)] #informative on betas

W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
#Special case: diagonal matrix
#(you may choose something else)
for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
    W[indexMoment,indexMoment] = 1.0/(dictEmpiricalMoments[k][1])^2
end

# Send to workers
sendto(workers(), dictPriors=dictPriors)
sendto(workers(), dictEmpiricalMoments=dictEmpiricalMoments)
sendto(workers(), W=W)

@everywhere set_priors!(myProblem, dictPriors)
@everywhere set_empirical_moments!(myProblem, dictEmpiricalMoments)
@everywhere set_weight_matrix!(myProblem, W)

# x[1] corresponds to the intercept, x[2] corresponds to beta1, x[3] corresponds to beta2
@everywhere function functionLinearModel(x; uniform_draws::Array{Float64,1}, simX::Array{Float64,2}, nbDraws::Int64 = length(uniform_draws), burnInPerc::Int64 = 0)
    T = nbDraws
    P = 2       #number of dependent variables

    alpha = x[1]
    beta = x[2:end]
    theta = 0.0     #coefficient to create serial correlation in the error terms

    # Creation of error terms
    # row = individual dimension
    # column = time dimension
    U = zeros(T)
    d = Normal()
    # Inverse cdf (i.e. quantile)
    gaussian_draws = quantile.(d, uniform_draws)
    U[1] = gaussian_draws[1] #first error term

    # loop over time periods
    for t = 2:T
        U[t] = gaussian_draws[t] + theta*U[t-1]
    end

    # Let's calculate the resulting y_t
    y = zeros(T)

    for t=1:T
        y[t] = alpha + simX[t,1]*beta[1] + simX[t,2]*beta[2] + U[t]
    end

    # Get rid of the burn-in phase:
    #------------------------------
    startT = max(1, Int(nbDraws * (burnInPerc / 100)))

    # Moments:
    #---------
    output = OrderedDict{String,Float64}()
    output["mean"] = mean(y[startT:nbDraws])
    output["mean^2"] = mean(y[startT:nbDraws].^2)
    output["mean_x1y"] = mean(simX[startT:nbDraws,1] .* y[startT:nbDraws])
    output["mean_x2y"] = mean(simX[startT:nbDraws,2] .* y[startT:nbDraws])
    output["mean_x1y^2"] = mean((simX[startT:nbDraws,1] .* y[startT:nbDraws]).^2)
    output["mean_x2y^2"] = mean((simX[startT:nbDraws,2] .* y[startT:nbDraws]).^2)

    return output
end

# Let's freeze the randomness during the minimization
d_Uni = Uniform(0,1)
nbDraws = T #Number of draws in the simulated data
burnInPerc = 0 #Burn-in phase (10%). Not necessary in the present context.
startT = max(1, Int(nbDraws * (burnInPerc / 100))) #First period used to calculate moments on simulated data
NMSM = nbDraws - startT + 1; #Number of Draws used when calculated moments on simulated data
uniform_draws = rand(d_Uni, nbDraws)
simX = zeros(length(uniform_draws), 2)
d = Uniform(0, 5)
for p = 1:2
  simX[:,p] = rand(d, length(uniform_draws))
end

# Send to workers
sendto(workers(), burnInPerc=burnInPerc)
sendto(workers(), simX=simX)
sendto(workers(), uniform_draws=uniform_draws)

# Construct the objective function everywhere
@everywhere set_simulate_empirical_moments!(myProblem, x -> functionLinearModel(x, uniform_draws = uniform_draws, simX = simX, burnInPerc=burnInPerc))
@everywhere construct_objective_function!(myProblem)

# Safety check: value on the master node == values on slave nodes?
using Test
val_local = myProblem.objective_function(ones(3)); #local execution
val_workers = [];
for w in workers() #Execution on workers
    push!(val_workers, @fetchfrom w myProblem.objective_function(ones(3)))
end
for (wIndex, w) in enumerate(workers())
    @test abs(val_local - val_workers[wIndex]) < 10e-10
end

#------------------------------------------------------------------------------
# Formulate the problem as a Laplace Type Estimator and use MCMC to find
# the quasi-posterior median
#------------------------------------------------------------------------------
# For tuning parameters, see: See https://github.com/madsjulia/AffineInvariantMCMC.jl
@everywhere begin
	numdims = 3
	numwalkers = 10
	thinning = 10
	numsamples_perwalker = 10000
	burnin = Int((1/10)*numsamples_perwalker)
	lb = 0 .* ones(numdims) #lower bound
	ub = 1 .* ones(numdims) #upper bound
	# Uniform prior
	# d_prior = Product(Uniform.(lb, ub))
	# Normal
	d_prior = MvNormal(zeros(numdims), 1.0 .* I(numdims))
end

# Pseudo Log-likelihood
@everywhere function Ln_MSM(x, NMSM)
	return -0.5*NMSM*myProblem.objective_function(x)
end

# Pseudo Log quasi-posterior: Pseudo Log(likelihood) + log(prior)
@everywhere function quasi_posterior(x, NMSM, d_prior)
	 return Ln_MSM(x, NMSM) + log(pdf(d_prior, x))
end

# Safety check: value on the master node == values on slave nodes?
using Test
val_local = Ln_MSM(ones(3), NMSM); #local execution
val_workers = [];
for w in workers() #Execution on workers
    push!(val_workers, @fetchfrom w Ln_MSM(ones(3), NMSM))
end
for (wIndex, w) in enumerate(workers())
    @test abs(val_local - val_workers[wIndex]) < 10e-10
end

# Slightly perturb the initial draws for the walkers
x0 = [dictPriors[k][1] for k in keys(dictPriors)]
x0_chains = ones(numdims, numwalkers).*true_vals .+ rand(numdims, numwalkers) .* 1.0
chain, llhoodvals = AffineInvariantMCMC.sample(x -> quasi_posterior(x, nbDraws, d_prior), numwalkers, x0_chains, burnin, 1)
chain, llhoodvals = AffineInvariantMCMC.sample(x -> quasi_posterior(x, nbDraws, d_prior), numwalkers, chain[:, :, end], numsamples_perwalker, thinning)
flatchain, flatllhoodvals = AffineInvariantMCMC.flattenmcmcarray(chain, llhoodvals)

#-------------------------------------------------------------------------------
# Plot Draws
#-------------------------------------------------------------------------------
p1 = plot(flatchain[1,:], ylabel="alpha0", xlabel="T", legend=:none)
p2 = plot(flatchain[2,:], ylabel="beta1", xlabel="T", legend=:none)
p3 = plot(flatchain[3,:], ylabel="beta2", xlabel="T", legend=:none)
p4 = plot(p1, p2, p3)
savefig(p4, joinpath(pwd(),"chains_MSM_MCMC.png"))

hh1 = histogram(flatchain[1,burnin:end], title="alpha0", legend=:none)
vline!(hh1, [alpha0[1]], linewidth = 4)
hh2 = histogram(flatchain[2,burnin:end], title="beta1", legend=:none)
vline!(hh2, [beta0[1]], linewidth = 4)
hh3 = histogram(flatchain[3,burnin:end], title="beta2", legend=:none)
vline!(hh3, [beta0[2]], linewidth = 4)
hh4 = plot(hh1, hh2, hh3)
savefig(hh4, joinpath(pwd(),"histograms_MSM_MCMC.png"))

# Compare results with GLM
using DataFrames, GLM
data = DataFrame(x1=x[:,1], x2=x[:,2], y= y[:]);
ols = lm(@formula(y ~ x1 + x2), data)
coef_ols = coef(ols)
ci_ols = confint(ols)
stderror_ols = stderror(ols)

result_alpha0 = append!(quantile(flatchain[1,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[1,burnin:end]), NaN, alpha0[1],NaN, coef_ols[1], ci_ols[1,1], ci_ols[1,2], stderror_ols[1])
result_beta1 = append!(quantile(flatchain[2,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[2,burnin:end]), NaN,beta0[1], NaN,coef_ols[2], ci_ols[2,1], ci_ols[2,2], stderror_ols[2])
result_beta2 = append!(quantile(flatchain[3,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[3,burnin:end]), NaN,beta0[2], NaN,coef_ols[3], ci_ols[3,1], ci_ols[3,2], stderror_ols[3])
results = DataFrame(variable = ["P5"; "P10"; "Median"; "P90"; "P95"; "std"; "-" ;"True value"; "-" ;"OLS Estimate"; "P5 OLS"; "P95 OLS"; "Std OLS"],
						alpha0 = result_alpha0, beta1 = result_beta1, beta2 = result_beta2)

CSV.write(joinpath(pwd(),"output_table_MSM_MCMC.csv"), results)

#-------------------------------------------------------------------------------
# Efficient MSM
#-------------------------------------------------------------------------------
# Set the W = inverse(Distance Matrix)
# Distance Matrix using Empirical Series
#---------------------------------------
X = zeros(T, length(dictEmpiricalMoments))
X[:,1] = y
X[:,2] = y.^2
X[:,3] = (x[:,1] .* y)
X[:,4] = (x[:,2] .* y)
X[:,5] = (x[:,1] .* y).^2
X[:,6] = (x[:,2] .* y).^2
Sigma0 = cov(X)
W_efficient = inv(Sigma0)
sendto(workers(), W_efficient=W_efficient)

@everywhere set_weight_matrix!(myProblem, W_efficient)
@everywhere construct_objective_function!(myProblem)

# Safety check: value on the master node == values on slave nodes?
using Test
val_local = myProblem.objective_function(ones(3)); #local execution
val_workers = [];
for w in workers() #Execution on workers
    push!(val_workers, @fetchfrom w myProblem.objective_function(ones(3)))
end
for (wIndex, w) in enumerate(workers())
    @test abs(val_local - val_workers[wIndex]) < 10e-10
end

#Slightly perturb the initial draws for the walkers
x0 = [dictPriors[k][1] for k in keys(dictPriors)]
x0_chains = ones(numdims, numwalkers).*true_vals .+ rand(numdims, numwalkers) .* 1.0
chain, llhoodvals = AffineInvariantMCMC.sample(x -> quasi_posterior(x, nbDraws, d_prior), numwalkers, x0_chains, burnin, 1)
chain, llhoodvals = AffineInvariantMCMC.sample(x -> quasi_posterior(x, nbDraws, d_prior), numwalkers, chain[:, :, end], numsamples_perwalker, thinning)
flatchain, flatllhoodvals = AffineInvariantMCMC.flattenmcmcarray(chain, llhoodvals)

#-------------------------------------------------------------------------------
# Plot Draws
#-------------------------------------------------------------------------------
p1 = plot(flatchain[1,:], ylabel="alpha0", xlabel="T", legend=:none)
p2 = plot(flatchain[2,:], ylabel="beta1", xlabel="T", legend=:none)
p3 = plot(flatchain[3,:], ylabel="beta2", xlabel="T", legend=:none)
p4 = plot(p1, p2, p3)
savefig(p4, joinpath(pwd(),"chains_MSM_MCMC_efficient.png"))

hh1 = histogram(flatchain[1,burnin:end], title="alpha0", legend=:none)
vline!(hh1, [alpha0[1]], linewidth = 4)
hh2 = histogram(flatchain[2,burnin:end], title="beta1", legend=:none)
vline!(hh2, [beta0[1]], linewidth = 4)
hh3 = histogram(flatchain[3,burnin:end], title="beta2", legend=:none)
vline!(hh3, [beta0[2]], linewidth = 4)
hh4 = plot(hh1, hh2, hh3)
savefig(hh4, joinpath(pwd(),"histograms_MSM_MCMC_efficient.png"))

# Compare results with GLM
using DataFrames, GLM
data = DataFrame(x1=x[:,1], x2=x[:,2], y= y[:]);
ols = lm(@formula(y ~ x1 + x2), data)
coef_ols = coef(ols)
ci_ols = confint(ols)
stderror_ols = stderror(ols)

result_alpha0 = append!(quantile(flatchain[1,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[1,burnin:end]), NaN, alpha0[1],NaN, coef_ols[1], ci_ols[1,1], ci_ols[1,2], stderror_ols[1])
result_beta1 = append!(quantile(flatchain[2,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[2,burnin:end]), NaN,beta0[1], NaN,coef_ols[2], ci_ols[2,1], ci_ols[2,2], stderror_ols[2])
result_beta2 = append!(quantile(flatchain[3,burnin:end],[0.05, 0.10, 0.5, 0.90, 0.95]), std(flatchain[3,burnin:end]), NaN,beta0[2], NaN,coef_ols[3], ci_ols[3,1], ci_ols[3,2], stderror_ols[3])
results_efficient = DataFrame(variable = ["P5"; "P10"; "Median"; "P90"; "P95"; "std"; "-" ;"True value"; "-" ;"OLS Estimate"; "P5 OLS"; "P95 OLS"; "Std OLS"],
					   alpha0 = result_alpha0, beta1 = result_beta1, beta2 = result_beta2)

CSV.write(joinpath(pwd(),"output_table_MSM_MCMC_efficient.csv"), results_efficient)


text:
### useful plotting functions

### savedisplay save svg and display current figure

function savefigure(label, fmts=[".svg", ".png"], projdir=getprojectdir(); p = plot!(), verbose=true)
    verbose && println("Saving Figure: ", label)
    files = joinpath.(projdir, "figures", string.(label, fmts))
    mkpath(dirname(first(files)))
    plot!(p)
    map(savefig, files)
end

function savedisplay(label, fmts=[".svg", ".png"], projdir=getprojectdir(); p = plot!(), verbose=false)
    savefigure(label, fmts, projdir, p=p, verbose=verbose)
    display(p)
end

######## mutation plots

function plotmutations(region, mutations)
    mutfam = mutations[!, Not([:Patient, :Inheritance])] |> unique
    mutgroups = combine(groupby(mutfam[mutfam.Contained, :], [:RelReg]), nrow => :Count, :type => Ref => :type, :Alt => Ref => :Alt)

    p = plotseq(region.metadata)

    acgt = MotifScanner.motif_letter_data(weight="normal")

    for row in eachrow(mutgroups)
        for i = 1:row.Count
            if row.type[i] == "var"
                plotletter!([row.RelReg[1] .- 0.5, -i*0.275], [1.0, 0.25], acgt[first(row.Alt[i])], c=:red) ### hard coded sizes to match MotifScanner, TODO: remove these
            end

            if row.type[i] == "del"
                plot!(row.RelReg, -0.25*ones(length(row.RelReg)), c=:red, lab="")
                plot!(row.RelReg[1]*[1, 1], [-0.1, -0.25], c=:red, lab="")
                plot!(row.RelReg[end]*[1, 1], [-0.1, -0.25], c=:red, lab="")
                annotate!([(mean(row.RelReg), -0.5, text("del", font("helvetica oblique", 9, :red)))])
            end
        end
    end
    p
end

"""
    idconvert(x)

    Function to trim lead 711086.. off snp Ids
"""
function idconvert(x)
    fields = split(x, "_")
    string(replace(fields[1], "711086" => ""), ":", fields[2], ">", fields[3])
end

function disruptmotifscatter(motsel ; kwargs...)

    theme(:wong2)

    @with motsel scatter(:RefPrMax, :AltPrMax, group=:MotifFamily, marker=(stroke(0)))#, marker=^(:auto))
    plot!([0, 1], [0, 1], c=:black, lab="", leg=:outertopright, xlabel="Motif Score Ref.", ylabel="Motif Score Mut")
    xlims!(0.5, 0.8)
    ylims!(-.125, 0.8)
    yticks!(-0.0:.2:0.8)
    annotate!([(motsel.RefPrMax[i] .+ 0.01, motsel.AltPrMax[i], text(motsel.ID[i] |> idconvert, "helvetica", :left, 6)) for i = 1:2:size(motsel, 1)])
    annotate!([(motsel.RefPrMax[i] .- 0.01, motsel.AltPrMax[i], text(motsel.ID[i] |> idconvert, "helvetica", :right, 6)) for i = 2:2:size(motsel, 1)])
    ps = plot!(size=(275, 225), grid=false, title="Top Ranked Motifs\ndisrupted by mutation", fontfamily="helvetica", titlefont=font(10, "helvetica"); kwargs...)
    theme(:wong)
    ps

end

function disruptedmotif_seqlogo(region, mutations, motsel)
    p = plotmutations(region, mutations)
    plot!(size=(600, 200), ylims=(-2, 3.5), grid=false, yaxis=false, yticks=false)

    for row in eachrow(motsel)
        motif = motifs[findfirst(motifmeta.MotifName .== row.MotifName)]
        xo = row.RefSeqStart-1
        yo = 0.5 #+ ifelse(row.MotifName .== "Hic1", 2, 0)
        seqlogo!(motif, rc=row.RefStrand == "-", xo=xo, yo=yo)
        annotate!((xo + size(motif.pwm, 2)/2, yo + 2.5, text(row.MotifFamily, font("helvetica", :top, :center, 8))))
    end
    xl = (71108635, 71108685+12)
    plot!(xticks=(1:25:101, string.( (1:25:101) .+ leftposition(region) .- 1)), xlims = xl .- leftposition(region) .+ 1, xlabel=string("chr10:", first(xl), "-", last(xl)))
    plot!(bottom_margin=5mm, fontfamily="helvetica")
    p
end

##### heatmap averaging

function avh(H::Vector{T}, δ) where {T}
    n = length(H)

    w = cld(n, δ)
    AH = zeros(Float64, w)
    te = zeros(Int, w)
    for i = 1:n
        wi = cld(i, δ)
        te[wi] += 1
        AH[wi] += H[i]
    end
    if te[end] < 0.75δ
        AH[end] += AH[end-1]
        te[end] += te[end-1]
    end
    AH./te
end

############### plotting gene models

rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
function plotgene!(trans, exoncoords;  y=1, off=1, h = 0.5, label="", fs=10, pos=:left, shiftstrand=true, c=:black, kwargs...)

    g = @subset(exoncoords, :TranscriptID .== trans)
    p = plot!(size=(1200, 100) ; kwargs...)

    if shiftstrand
        if g.strand[1] == "-"
            # y = y - 1
            # c = :red
        end
    end

    for (s, e) in g.coords[1]
       plot!(rectangle(e - s + 1, h, s - off + 1, y), c=c, linecolor=c, lab="")
    end
    plot!([g.start[1], g.stop[1]] .- off .+ 1, y .+ h/2 .+ [0, 0], c=c, lab="")
    if !isempty(label)
        if pos == :left
            annotate!([(g.start[1] - 2000, y + h/2, text(label, font(:right, fs, "helvetica")))])
        else
            ### corrections to yoff
            yc = 0
            if g.GeneName[1] == "HK1"
                yc = -1
            elseif g.GeneName[1] .== "TYSND1"
                yc = 1.5
            end
            annotate!([((g.start[1] + g.stop[1])/2, y-0.5 + yc, text(label,font(:center, fs, "helvetica")))])
        end
    end
    p
end

overlaps(startA, stopA, interval) = !isempty(intersect(startA:stopA, interval))

locs(s, e) = s:e
function plotgenemodels(loc, exoncoords; off=1, fs=10)

    maxmodels = @chain exoncoords begin
        @subset(overlaps.(:start, :stop, Ref(loc)))
        groupby(:GeneName)
        combine(df -> df[argmax(df.stop - df.start), :])
    end

    maxmodels.ovg = overlappinglocations(fill("", size(maxmodels, 1)), locs.(maxmodels.start .- 1000, maxmodels.stop .+ 1000))
    maxmodels = transform(groupby(maxmodels, :ovg), :ovg => (x -> 1:length(x)) => :yoff)
    maxmodels[maxmodels.GeneName .== "TACR2", :yoff] .+= 1
    p  = plot()

    for (g, t, yo) in zip(maxmodels.GeneName, maxmodels.TranscriptID, maxmodels.yoff)
        plotgene!(t, exoncoords, y=1 - yo + 1, off=off, label=g, fs=fs, pos=:bottom)
    end
    if "ENST00000359426.7_2" ∈ exoncoords.TranscriptID
        plotgene!("ENST00000359426.7_2", exoncoords, y=0, off=off, label="", fs=fs, pos=:bottom) ## extra isoform for hk1 shortest and highested expressed in islets
        #plotgene!("ENST00000298649.8_2", exoncoords, y=-1, off=off, label="", fs=fs, pos=:bottom) ##  second shortest, there is TF binding right on the promoter of this, but does not appear to be accessible or expressed
    end

    xlims!(first(loc), last(loc))
    p
end


text:
# Pkg.test runs with --check_bounds=1, forcing all bounds checks.
# This is incompatible with GPUifyLoops.
if Base.JLOptions().check_bounds == 1
  cmd = Cmd(filter(arg->!startswith(arg, "--check-bounds"), Base.julia_cmd().exec))
  code = """
    $(Base.load_path_setup_code(false))
    cd($(repr(@__DIR__)))
    include($(repr(@__FILE__)))
    """
  run(`$cmd --eval $code`)
  exit()
end
@assert Base.JLOptions().check_bounds == 0

using SafeTestsets, Test

@time @safetestset "EnsembleGPUArray" begin include("ensemblegpuarray.jl") end
@time @safetestset "EnsembleGPUArray SDE" begin include("ensemblegpuarray_sde.jl") end
@time @safetestset "EnsembleGPUArray Input Types" begin include("ensemblegpuarray_inputtypes.jl") end
# Not safe because distributed doesn't play nicely with modules.
@time @testset "Distributed Multi-GPU" begin include("distributed_multi_gpu.jl") end


text:
using Clang.Generators
using Clingo_jll

cd(@__DIR__)

const HEADER_BASE = joinpath(Clingo_jll.artifact_dir, "include")
const CLINGO_H = joinpath(HEADER_BASE, "clingo.h")

headers = [CLINGO_H]

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = ["-I$HEADER_BASE"]

ctx = create_context(headers, args, options)

build!(ctx)


text:
using Documenter
using ReferenceFrameRotations

DocMeta.setdocmeta!(
    ReferenceFrameRotations,
    :DocTestSetup,
    :(using ReferenceFrameRotations);
    recursive=true
)

makedocs(
    modules = [ReferenceFrameRotations],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://juliaspace.github.io/ReferenceFrameRotations.jl/stable/",
    ),
    sitename = "Reference Frame Rotations",
    authors = "Ronan Arraes Jardim Chagas",
    pages = [
        "Home" => "index.md",
        "Direction Cosine Matrices" => "man/dcm.md",
        "Euler Angle and Axis" => "man/euler_angle_axis.md",
        "Euler Angles" => "man/euler_angles.md",
        "Quaternions" => "man/quaternions.md",
        "Conversions" => "man/conversions.md",
        "Kinematics" => "man/kinematics.md",
        "Composing rotations" => "man/composing_rotations.md",
        "Inverting rotations" => "man/inv_rotations.md",
        "Library" => "lib/library.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaSpace/ReferenceFrameRotations.jl.git",
    target = "build",
)


text:
include("im_io/matlab.jl")
include("im_io/common.jl")
include("im_io/data.jl")


text:
# - Prerequesites:
# --- either:
# ----- download wahlomat data manually from:
#       https://www.bpb.de/system/files/datei/Wahl-O-Mat%20Bundestag%202021_Datensatz_v1.02.zip
# ----- unzip, rename to 'wahlomat.xlsx', delete everything else
# --- or:
# ----- ./get_data.sh

# Activate julia env
using Pkg
Pkg.activate(".")
Pkg.update()

# Reset database
if isfile("./db.sqlite")
    Base.Filesystem.rm("./db.sqlite")  # hard reset: if db exists, drop it
end

# Dependencies for formatting
using Negotiations
using DataFrames
using XLSX
using SQLite
using Match

# function to code opinion numerically
function recode_opinion(opinion::String)
    @match opinion begin
        "stimme zu" => 1
        "neutral" => 0
        "stimme nicht zu" => -1
        _ => -10
    end
end

printstyled("  Processing "; color = :green, bold = true)
print("data...\n")

# Read data
df = DataFrame(XLSX.readtable("data/wahlomat.xlsx", "Datensatz Bundestag 2021")...)

# Sane column names
rename!(
    df,
    [
        "party_id",
        "party_shorthand",
        "party_name",
        "statement_id",
        "statement_title",
        "statement",
        "position",
        "position_rationale"
    ]
)

# Split and process tables
party_table = unique(df[:, [:party_id, :party_shorthand, :party_name]])
party_table.party_shorthand = replace.(party_table.party_shorthand, "Ä" => "AE")
party_table.party_shorthand = replace.(party_table.party_shorthand, "ä" => "ae")
party_table.party_shorthand = replace.(party_table.party_shorthand, "Ö" => "OE")
party_table.party_shorthand = replace.(party_table.party_shorthand, "ö" => "oe")
party_table.party_shorthand = replace.(party_table.party_shorthand, "Ü" => "UE")
party_table.party_shorthand = replace.(party_table.party_shorthand, "ü" => "ue")
party_table.party_shorthand = replace.(party_table.party_shorthand, "/" => " ")
party_table.party_shorthand = replace.(party_table.party_shorthand, "³" => "3")
party_table.party_shorthand = replace.(party_table.party_shorthand, r"\." => " ")
party_table.party_shorthand = replace.(party_table.party_shorthand, r"-" => " ")
party_table.party_shorthand = replace.(party_table.party_shorthand, " " => "")

statement_table = unique(df[:, [:statement_id, :statement_title, :statement]])

opinion_table = df[:, [:party_id, :statement_id, :position, :position_rationale]]
opinion_table.position = [recode_opinion(o) for o in opinion_table.position]

printstyled("  Creating "; color = :green, bold = true)
print("database...\n")

# Create database
db = initialize_db("db")
SQLite.load!(party_table, db, "party")
SQLite.load!(statement_table, db, "statement")
SQLite.load!(opinion_table, db, "opinion")


text:
using NeXLParticle
using StaticArrays
using GeometryBasics
using Random
using LinearAlgebra
using Rotations
using Gadfly

r0 = MersenneTwister(0xEA7F00D)
for npts in (100, 1000, 10000, 100000)
    println("Working on $npts points...")
    for i in 1:10
        seed = rand(r0, UInt)
        rnd = MersenneTwister(seed)
        snp = sqrt(npts)
        bounds = Rect2((-0.5*snp, -0.5*snp), (snp,snp))
        xy = NeXLParticle.generate_ground_truth(rnd, bounds, npts, x->norm(x)<0.5*snp)
        θ1, θ2 = 2π*rand(rnd, 2)
        f1, f2 = 0.5 .+ 0.5 * rand(rnd, 2)
        ps1 = NeXLParticle.measure_particles(rnd, xy, SA[100.0*rand(rnd),100.0*rand(rnd)], θ1, f1, 0.001)
        ps2 = NeXLParticle.measure_particles(rnd, xy, SA[100.0*rand(rnd),100.0*rand(rnd)], θ2, f2, 0.001)

        (ct1, ct2) = @time rough_align(ps1, ps2, tol=0.001)

        r = RotMatrix{2}(θ1-θ2)
        if !isapprox(r, ct2.linear, atol = 0.01)

            7println("* * * * rough_align(...) failure * * * *")
            println(" Failure:  $seed")
            println(" NPoints: $npts")
            println("Fraction:f1 = $f1, f2 = $f2")
            println("   Angle: θ1 = $θ1, θ2 = $θ2")
            println("   Truth: $r")
            println("Estimate: $ct2.linear")
            ts1, ts2 = ct1.(ps1), ct2.(ps2)
            p = plot( ( layer(x=map(p->p[1], t), y=map(p->p[2], t), Geom.point, Theme(default_color=c)) for (t,c) in zip((ts1, ts2), ("red", "blue"))) ...)
            p |> SVG(joinpath(homedir(),"Desktop","RA[$seed].svg"), 20inch, 20inch)
        end
        GC.gc()
    end
end
#=
The algorithm seems to fail when:
  1) the number of datapoints is small and the measurement fraction is low
  2) density of particles is too high
  3) variability in measured positions is high

  Offset and angles don't matter.

Timings:
Working on 100 points... (Average: 4.3 ms)
  0.002455 seconds (34.60 k allocations: 1.318 MiB)
  0.002719 seconds (56.89 k allocations: 2.134 MiB)
  0.003065 seconds (63.53 k allocations: 2.380 MiB)
  0.001486 seconds (29.57 k allocations: 1.121 MiB)
  0.006415 seconds (134.60 k allocations: 5.160 MiB)
  0.004236 seconds (67.86 k allocations: 2.573 MiB)
  0.004809 seconds (104.54 k allocations: 3.893 MiB)
  0.004412 seconds (51.79 k allocations: 1.950 MiB)
  0.006890 seconds (139.50 k allocations: 5.320 MiB)
  0.006069 seconds (113.35 k allocations: 4.214 MiB)
Working on 1,000 points...  (Average: 9.99 ms)
  0.010888 seconds (185.80 k allocations: 7.901 MiB)
  0.013183 seconds (183.53 k allocations: 7.804 MiB)
  0.011232 seconds (182.69 k allocations: 7.709 MiB)
  0.008922 seconds (180.93 k allocations: 7.608 MiB)
  0.010046 seconds (182.18 k allocations: 7.679 MiB)
  0.006969 seconds (134.49 k allocations: 5.283 MiB)
  0.010058 seconds (185.00 k allocations: 7.854 MiB)
  0.008368 seconds (158.51 k allocations: 6.278 MiB)
  0.010471 seconds (185.53 k allocations: 7.887 MiB)
  0.009765 seconds (181.95 k allocations: 7.672 MiB)
Working on 10,000 points...  (Average: 29.5 ms)
  0.037703 seconds (360.80 k allocations: 18.659 MiB)
  0.027629 seconds (305.02 k allocations: 15.226 MiB)
  0.027257 seconds (291.41 k allocations: 14.402 MiB)
  0.027362 seconds (310.73 k allocations: 15.577 MiB)
  0.025490 seconds (300.53 k allocations: 14.935 MiB)
  0.026743 seconds (297.59 k allocations: 14.743 MiB)
  0.034434 seconds (353.36 k allocations: 18.208 MiB)
  0.034588 seconds (354.96 k allocations: 18.272 MiB)
  0.031228 seconds (320.44 k allocations: 16.160 MiB)
  0.022824 seconds (269.03 k allocations: 13.118 MiB)
Working on 100,000 points...  (Average:  285.3 ms)
  0.270818 seconds (1.52 M allocations: 89.571 MiB, 6.31% gc time)
  0.349587 seconds (1.98 M allocations: 118.098 MiB, 5.10% gc time)
  0.231989 seconds (1.49 M allocations: 87.577 MiB)
  0.335018 seconds (1.89 M allocations: 112.289 MiB, 5.34% gc time)
  0.256628 seconds (1.50 M allocations: 88.605 MiB, 7.25% gc time)
  0.283030 seconds (1.64 M allocations: 96.868 MiB, 6.65% gc time)
  0.298009 seconds (1.64 M allocations: 96.943 MiB, 7.23% gc time)
  0.348799 seconds (2.01 M allocations: 120.019 MiB, 5.38% gc time)
  0.221129 seconds (1.37 M allocations: 80.252 MiB)
  0.257807 seconds (1.51 M allocations: 89.376 MiB, 6.76% gc time)
=#

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: sagemaker_a2i_runtime
using AWS.Compat
using AWS.UUIDs

"""
    delete_human_loop(human_loop_name)
    delete_human_loop(human_loop_name, params::Dict{String,<:Any})

Deletes the specified human loop for a flow definition. If the human loop was deleted, this
operation will return a ResourceNotFoundException.

# Arguments
- `human_loop_name`: The name of the human loop that you want to delete.

"""
delete_human_loop(HumanLoopName; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("DELETE", "/human-loops/$(HumanLoopName)"; aws_config=aws_config)
delete_human_loop(HumanLoopName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("DELETE", "/human-loops/$(HumanLoopName)", params; aws_config=aws_config)

"""
    describe_human_loop(human_loop_name)
    describe_human_loop(human_loop_name, params::Dict{String,<:Any})

Returns information about the specified human loop. If the human loop was deleted, this
operation will return a ResourceNotFoundException error.

# Arguments
- `human_loop_name`: The name of the human loop that you want information about.

"""
describe_human_loop(HumanLoopName; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("GET", "/human-loops/$(HumanLoopName)"; aws_config=aws_config)
describe_human_loop(HumanLoopName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("GET", "/human-loops/$(HumanLoopName)", params; aws_config=aws_config)

"""
    list_human_loops(flow_definition_arn)
    list_human_loops(flow_definition_arn, params::Dict{String,<:Any})

Returns information about human loops, given the specified parameters. If a human loop was
deleted, it will not be included.

# Arguments
- `flow_definition_arn`: The Amazon Resource Name (ARN) of a flow definition.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"CreationTimeAfter"`: (Optional) The timestamp of the date when you want the human loops
  to begin in ISO 8601 format. For example, 2020-02-24.
- `"CreationTimeBefore"`: (Optional) The timestamp of the date before which you want the
  human loops to begin in ISO 8601 format. For example, 2020-02-24.
- `"MaxResults"`: The total number of items to return. If the total number of available
  items is more than the value specified in MaxResults, then a NextToken is returned in the
  output. You can use this token to display the next page of results.
- `"NextToken"`: A token to display the next page of results.
- `"SortOrder"`: Optional. The order for displaying results. Valid values: Ascending and
  Descending.
"""
list_human_loops(FlowDefinitionArn; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("GET", "/human-loops", Dict{String, Any}("FlowDefinitionArn"=>FlowDefinitionArn); aws_config=aws_config)
list_human_loops(FlowDefinitionArn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("GET", "/human-loops", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("FlowDefinitionArn"=>FlowDefinitionArn), params)); aws_config=aws_config)

"""
    start_human_loop(flow_definition_arn, human_loop_input, human_loop_name)
    start_human_loop(flow_definition_arn, human_loop_input, human_loop_name, params::Dict{String,<:Any})

Starts a human loop, provided that at least one activation condition is met.

# Arguments
- `flow_definition_arn`: The Amazon Resource Name (ARN) of the flow definition associated
  with this human loop.
- `human_loop_input`: An object that contains information about the human loop.
- `human_loop_name`: The name of the human loop.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"DataAttributes"`: Attributes of the specified data. Use DataAttributes to specify if
  your data is free of personally identifiable information and/or free of adult content.
"""
start_human_loop(FlowDefinitionArn, HumanLoopInput, HumanLoopName; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("POST", "/human-loops", Dict{String, Any}("FlowDefinitionArn"=>FlowDefinitionArn, "HumanLoopInput"=>HumanLoopInput, "HumanLoopName"=>HumanLoopName); aws_config=aws_config)
start_human_loop(FlowDefinitionArn, HumanLoopInput, HumanLoopName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("POST", "/human-loops", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("FlowDefinitionArn"=>FlowDefinitionArn, "HumanLoopInput"=>HumanLoopInput, "HumanLoopName"=>HumanLoopName), params)); aws_config=aws_config)

"""
    stop_human_loop(human_loop_name)
    stop_human_loop(human_loop_name, params::Dict{String,<:Any})

Stops the specified human loop.

# Arguments
- `human_loop_name`: The name of the human loop that you want to stop.

"""
stop_human_loop(HumanLoopName; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("POST", "/human-loops/stop", Dict{String, Any}("HumanLoopName"=>HumanLoopName); aws_config=aws_config)
stop_human_loop(HumanLoopName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = sagemaker_a2i_runtime("POST", "/human-loops/stop", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("HumanLoopName"=>HumanLoopName), params)); aws_config=aws_config)


text:
module Test15

import Automa
import Automa.RegExp: @re_str
using Test

@testset "Test15" begin
    a = re"a+"
    a.actions[:enter] = [:enter]
    a.actions[:all]   = [:all]
    a.actions[:final] = [:final]
    a.actions[:exit]  = [:exit]
    b = re"b+"
    b.actions[:enter] = [:enter]
    b.actions[:all]   = [:all]
    b.actions[:final] = [:final]
    b.actions[:exit]  = [:exit]
    ab = Automa.RegExp.cat(a, b)

    machine = Automa.compile(ab)
    last, actions = Automa.execute(machine, "ab")
    @test last == 0
    @test actions == [:enter, :all, :final, :exit, :enter, :all, :final, :exit]

    for generator in (:table, :inline, :goto), checkbounds in (true, false), clean in (true, false)
        ctx = Automa.CodeGenContext(generator=generator, checkbounds=checkbounds, clean=clean)
        validate = @eval function (data)
            logger = Symbol[]
            $(Automa.generate_init_code(ctx, machine))
            p_end = p_eof = sizeof(data)
            $(Automa.generate_exec_code(ctx, machine, :debug))
            return logger, cs == 0 ? :ok : cs < 0 ? :error : :incomplete
        end
        @test validate(b"ab") == ([:enter, :all, :final, :exit, :enter, :all, :final, :exit], :ok)
    end
end

end


text:
using Test, SafeTestsets

println("importing QBase.jl")
@time using QBase

println("Testing QBase.jl")
@time @testset "./src/QBase.jl"  begin

println("Testing Math Utilities : ")
@time @testset "Math Utilities" begin
    println("./test/unit/math/matrices.jl")
    @time @safetestset "./test/unit/math/matrices.jl" begin include("unit/math/matrices.jl") end

    println("./test/unit/math/combinatorics.jl")
    @time @safetestset "./test/unit/math/combinatorics.jl" begin include("unit/math/combinatorics.jl") end
end

println("Testing Types : ")
@time @testset "Types" begin
    println("./test/unit/types/probabilities.jl")
    @time @safetestset "./test/unit/types/probabilities.jl" begin include("unit/types/probabilities.jl") end

    println("./test/unit/types/brakets.jl")
    @time @safetestset "./test/unit/types/brakets.jl" begin include("unit/types/brakets.jl") end

    println("./test/unit/types/states.jl")
    @time @safetestset "./test/unit/types/states.jl" begin include("unit/types/states.jl") end

    println("./test/unit/types/unitaries.jl")
    @time @safetestset "./test/unit/types/unitaries.jl" begin include("unit/types/unitaries.jl") end

    println("./test/unit/types/measurements.jl")
    @time @safetestset "./test/unit/types/measurements.jl" begin include("unit/types/measurements.jl") end
end

println("Testing Constructors : ")
@time @testset "Constructors" begin
    println("./test/unit/constructors/brakets.jl")
    @time @safetestset "./test/unit/constructors/brakets.jl" begin include("unit/constructors/brakets.jl") end

    println("./test/unit/constructors/states.jl")
    @time @safetestset "./test/unit/constructors/states.jl" begin include("unit/constructors/states.jl") end

    println("./test/unit/constructors/unitaries.jl")
    @time @safetestset "./test/unit/constructors/unitaries.jl" begin include("unit/constructors/unitaries.jl") end

    println("./test/unit/constructors/measurements.jl")
    @time @safetestset "./test/unit/constructors/measurements.jl" begin include("unit/constructors/measurements.jl") end
end

println("Testing Main : ")
@time @testset "Main" begin
    println("./test/unit/evolve.jl")
    @time @safetestset "./test/unit/evolve.jl" begin include("unit/evolve.jl") end

    println("./test/unit/measure.jl")
    @time @safetestset "./test/unit/measure.jl" begin include("unit/measure.jl") end

    println("./test/unit/channels.jl")
    @time @safetestset "./test/unit/channels.jl" begin include("unit/channels.jl") end

    println("./test/unit/information.jl")
    @time @safetestset "./test/unit/information.jl" begin include("unit/information.jl") end
end

end


text:
function lowercase_drive(a)
    if length(a) >= 2 && a[2] == ':'
        return lowercase(a[1]) * a[2:end]
    else
        return a
    end
end

const SRC_DIR = joinpath(Sys.BINDIR, "..", "..", "base")
const RELEASE_DIR = joinpath(Sys.BINDIR, "..", "share", "julia", "base")
basepath(file) = normpath(joinpath((@static isdir(SRC_DIR) ? SRC_DIR : RELEASE_DIR), file))

maybe_quote(x) = (isa(x, Expr) || isa(x, Symbol)) ? QuoteNode(x) : x

# all calls in the current frame and the same "line"
calls_on_line(state, args...) = calls_on_line(state.frame, args...)
calls_on_line(::Nothing, args...) = []
function calls_on_line(frame::JuliaInterpreter.Frame, line=nothing)
    if line === nothing
        loc = JuliaInterpreter.whereis(frame)
        if loc === nothing
            return []
        end
        line = loc[2]
    end

    src = frame.framecode.src

    exprs = []
    for pc in frame.pc:length(src.codelocs)
        loc = JuliaInterpreter.whereis(frame, pc)

        if loc === nothing || loc[2] > line
            return exprs
        end

        expr = deepcopy(JuliaInterpreter.pc_expr(src, pc))
        if Meta.isexpr(expr, :call)
            for i in 1:length(expr.args)
                expr.args[i] = maybe_quote(expr.args[i])
            end
            push!(exprs, (pc = pc, expr = string("%", pc, " = ", prettyprint_expr(expr, src))))
        elseif Meta.isexpr(expr, :(=))
            expr = expr.args[2]
            push!(exprs, (pc = pc, expr = prettyprint_expr(expr, src)))
        end
    end
    exprs
end

function prettyprint_expr(expr, src)
    io = IOBuffer()
    prettyprint_expr(io, expr, src)
    return String(take!(io))
end

function prettyprint_expr(io, expr, src)
    if Meta.isexpr(expr, :call)
        for (i, ex) in enumerate(expr.args)
            if ex isa QuoteNode
                print(io, ex.value)
            elseif ex isa JuliaInterpreter.SlotNumber || ex isa Core.SlotNumber
                print(io, src.slotnames[ex.id])
            elseif ex isa JuliaInterpreter.SSAValue || ex isa Core.SSAValue
                print(io, '%', ex.id)
            else
                prettyprint_expr(io, ex, src)
            end
            if i == 1
                print(io, '(')
            elseif i == length(expr.args)
                print(io, ')')
            else
                print(io, ", ")
            end
        end
    else
        show(io, expr)
    end
end

# all calls in the current frame after (and including) the current pc
calls_in_frame(state) = calls_in_frame(state.frame)
calls_in_frame(::Nothing) = []
function calls_in_frame(frame::JuliaInterpreter.Frame)
    exprs = []
    for pc in frame.pc:length(frame.framecode.src.codelocs)
        expr = JuliaInterpreter.pc_expr(frame, pc)
        if Meta.isexpr(expr, :call)
            push!(exprs, (pc = pc, expr = expr))
        elseif Meta.isexpr(expr, :(=))
            expr = expr.args[2]
            push!(exprs, (pc = pc, expr = expr))
        end
    end
    exprs
end


text:
using AmbientForcing
using PowerDynamics, Distributions

# Lets create a small test power grid with three nodes
line_list = []
append!(line_list, [StaticLine(from = 1, to = 2, Y = -1im / 0.02 + 4)])
append!(line_list, [StaticLine(from = 1, to = 3, Y = -1im / 0.02 + 4)])

# The PQAlgebraic Node is our constraint. The Power Output of node 3 is fixed.
node_list = []
append!(node_list, [SlackAlgebraic(U = 1, Y_n = 0)])
append!(node_list, [FourthOrderEq(H = 3.318, P = -0.6337, D = 0.1, Ω = 50, E_f = 0.5, T_d_dash = 0.1,  T_q_dash = 8.690,X_d_dash = 0.111,  X_q_dash = 0.103, X_d = 0.1, X_q = 0.6)])
append!(node_list, [PQAlgebraic(P = -0.6337, Q = 0.0)])

# Using the rigth hand side as our ODEFunction
pg = PowerGrid(node_list, line_list)
rpg = rhs(pg)

# Acessing the constraint equation g of the power grid
g = constraint_equations(rpg)

# The operation point is a fixed point and naturally lies on the manifold
# Thus it can be used as the initial condition for our differential equation
op = find_operationpoint(pg)

# Lets check if g(op) ≈ 0 meaning that the constraint is fulfilled
g(op.vec)

# Lets generate a random vector from the ambient space
# First we want to perturb all variables in the grid

Frand = random_force(rpg, [0.0,1], Uniform)

# As we can see our new valid inital condition fulfills the constraints :-)
z_new = ambient_forcing(rpg, op.vec, 2.0, Frand)
g(z_new)

## Next: let's perturb just the voltage at node 2!
# Getting the index of the real and imaginary part of the voltage
idx = idx_exclusive(rpg, ["u_r_2", "u_i_2"])

# Now we generate a vector with only non-vanishing componants at the voltage of node 2
Frand = random_force(rpg, [0.0, 1.0], Uniform, idx)

# We solve the differential equation and find that the constraint remains fulfilled!
z_new = ambient_forcing(rpg, op.vec, 2.0, Frand);
g(z_new)

# It is also possible to perturb the variable using different distributions:
# Typically SNBS the angle θ and ω are pertubed differently
idx = idx_exclusive(rpg, ["u_r_2", "u_i_2", "θ_2", "ω_2"])
τ = 2.0 # the integration time

# When you want to sample eg. the angle θ from a box of [0, 2π]
# You have to make sure to devide the distribution argument by the integration time τ
dist_vec = [[0,1], [0, 2], [0, 2π] ./τ, [-5, 5] ./τ]

# Then everything works as beforehand. And we see again that we have created a valid initial condition!
Frand = Frand = random_force(rpg, dist_vec, Uniform, idx)
z_new = ambient_forcing(rpg, op.vec, τ, Frand)
g(z_new)


text:
# Used for adjusting strength of pixel.
function g(x, max_rgb)
    return 1 - (x / max_rgb)
end

# find maximum colour vector
function find_maximum(rgb::AbstractArray)
    max_rgb = 0
    for i in CartesianIndices(rgb)
        current_RGB = sqrt(rgb[i].r^2 + rgb[i].g^2  + rgb[i].b^2)
        if current_RGB > max_rgb
            max_rgb = current_RGB
        end
    end
    return max_rgb
end

function count_enemy(labels, target_pixel::CartesianIndex, N₈)
    enemy_count = 0
    @inbounds for j in CartesianIndices(N₈)
         if j != target_pixel && labels[j] != labels[target_pixel] && labels[j] != 0
            enemy_count += 1
        end
    end
    return enemy_count
end

# default rule for GrowCut() uses decreasing function and relative difference in
# colour to adjust strengths.
function colour_differance(θ::AbstractArray, rgb::AbstractArray, target_pixel::CartesianIndex, attacking_pixel::CartesianIndex, extras)
    Cdiff = rgb[target_pixel] - rgb[attacking_pixel]
    absCdiff = sqrt(Cdiff.r^2 + Cdiff.g^2 + Cdiff.b^2)
    adjusted_strength = g(absCdiff, extras) * θ[attacking_pixel]
    return adjusted_strength
end

# std update of pixel
function update_pixel(θ::AbstractArray, θₜ₊₁::AbstractArray, labels::AbstractArray, labelsₜ₊₁::AbstractArray, feature_vector::AbstractArray, E::AbstractArray, target_pixel::CartesianIndex, extras, t1, count, N₈, changedₜ₊₁, calculate_adjusted_strength::Function)
    @inbounds for j in CartesianIndices(N₈)
        adjusted_strength = calculate_adjusted_strength(θ, feature_vector, target_pixel, j, extras)
        change = false
        if adjusted_strength > θ[target_pixel] && E[j] < t1 && j != target_pixel
            labelsₜ₊₁[target_pixel] = labels[j]
            θₜ₊₁[target_pixel] = adjusted_strength
            count += 1
            if change == false
                push!(changedₜ₊₁, N₈)
                change = true
            end
        end
    end
    return θₜ₊₁, labelsₜ₊₁, count, changedₜ₊₁
end

# update of pixel if surrounded by too many enemies
function occupy_pixel(θ, θₜ₊₁, labels, labelsₜ₊₁, feature_vector, target_pixel::CartesianIndex, extras, count, N₈, changedₜ₊₁, calculate_adjusted_strength)
    θmax = 1.0
    location = target_pixel
    ltemp = labels[target_pixel]
    @inbounds for j in CartesianIndices(N₈)
        if j != target_pixel && θ[j] < θmax && labels[j] != labels[target_pixel]
            θmax = θ[j]
            location = j
        end
    end
    θₜ₊₁[target_pixel] = calculate_adjusted_strength(θ, feature_vector, target_pixel, location, extras)
    labelsₜ₊₁[target_pixel] = labels[location]
    count += 1
    push!(changedₜ₊₁, N₈)
    return θₜ₊₁, labelsₜ₊₁, count, changedₜ₊₁
end

# sets initial strengths
function set_strength(label)
    θ = zeros(Float64, axes(label))
    for i in CartesianIndices(label)
        if label[i] > 0
            θ[i] = 1.0
        end
    end
    return θ
end

"""
```
help_rule_template()
```

Outputs the following template for a new rule into the REPL:

```julia
function your_rule(strengths::Array{Float64,1}, [your feature_vector]::AbstractArray, target_pixel::CartesianIndex, attacking_pixel::CartesianIndex, [your additional_values]::Any)
   # your rule goes here
   return adjusted_strength
end

function your_extra_stuff([your feature_vector]::AbstractArray
   #your calculations for any additional values goes here
   #if you need more then one value they can be in an array or wrapped into a tuple
   return additional_values
end
```

"""
function help_rule_template()
    println("function your_rule(strengths::Array{Float64,1}, [your feature_vector]::AbstractArray, target_pixel::CartesianIndex, attacking_pixel::CartesianIndex, [your additional_values]::Any)")
    println("   #your rule goes here")
    println("   return adjusted_strength")
    println("end")
    println("")
    println("function your_extra_stuff([your feature_vector]::AbstractArray")
    println("   #your calculations for any additional values goes here")
    println("   #if you need more then one value they can be in an array or wrapped into a tuple")
    println("   return additional_values")
    println("end")
    println("")
    println("Use create_rule() to print to new file")
end

"""
```
create_rule(file::AbstractString="newrule.jl")
```

Outputs template for a new rule into a file.

# Arguments

The function arguments are described in more detail below.

##  `name_of_rule`
An `AbstractString` specifying the file to put the function in. Note: this can
override file contents, recomended to create a new file.

"""
function create_rule(file::AbstractString="newrule.jl")
    write(file,"function your_rule(strengths::Array{Float64,1}, [your feature_vector]::AbstractArray, target_pixel::CartesianIndex, attacking_pixel::CartesianIndex, [your additional_values]::Any)\n",
    "   #your rule goes here\n",
    "   return adjusted_strength\n",
    "end\n",
    "\n",
    "function your_extra_stuff([your feature_vector]::AbstractArray\n",
    "   #your calculations for any additional values goes here\n",
    "   #if you need more then one value they can be in an array or wrapped into a tuple\n",
    "   return additional_values\n",
    "end")
end


text:
type Vertex <: Container
	# Represents a Vertex within a Property Graph

	id::UUID.Uuid
	containedobject::Any
	attachedproperties::Dict{String,Any}

	graph::Any
	# incomingedges and outgoingEdgens are declared Set{Container} instead of Set{Edge}
	# in order to avoid mutually circular type declaration with Edge
	# see: https://github.com/JuliaLang/julia/issues/269
	incomingedges::Set{Container}
	outgoingedges::Set{Container}
	typelabel::String
	partiallyloaded::Bool

	function Vertex(properties::Dict{String,Any}=Dict{String,Any}())
		# Construct a Vertex with a set of initial property values

		v = new()
		v.id = UUID.v4()

		v.attachedproperties = properties
		v.containedobject = UnspecifiedValue
		v.typelabel = "Unspecified"

		v.incomingedges = Set{Container}()
		v.outgoingedges = Set{Container}()

		v.partiallyloaded = false

		return v
	end

	function Vertex(typelabel::String, properties::Dict{String,Any})
		# Construct a Vertex with a typelabel and a set of initial property values

		v = Vertex(properties)
		v.typelabel = typelabel

		return v
	end

	function Vertex(typelabel::String)
		# Construct a Vertex with a typelabel

		v = Vertex(typelabel, Dict{String,Any}())

		return v
	end
end

function vertexforobject(object::Any, properties::Dict{String,Any}=Dict{String,Any}())
	# Construct a Vertex for an existing object, and with a set of initial property values

	v = Vertex(properties)
	v.containedobject = object

	return v
end

function vertexforobject(typelabel::String, object::Any, properties::Dict{String,Any}=Dict{String,Any}())
	# Construct a Vertex with a type label, a reference to an existing object, and a set of initial property values

	v = Vertex(properties)
	v.typelabel = typelabel
	v.containedobject = object

	return v
end


text:
module SerializedElementArrays

using Serialization

# tempname only supports cleanup since Julia 1.4
function cleanup_default()
  return VERSION >= v"1.4"
end

# Like `tempname()` but compatible with Julia versions prior to 1.4
function temppath(path::AbstractString=tempdir(); cleanup=cleanup_default())
  if VERSION < v"1.4"
    if cleanup
      error(
        "Previous to Julia 1.4, `tempname` doesn't support the `cleanup` keyword argument so `cleanup` must be set to false.",
      )
    end
    if path ≠ tempdir()
      error(
        "`path` was specified as $path. Previous to Julia 1.4, `tempname` doesn't support customizing the path, `path` must be `tempdir() = $(tempdir())`.",
      )
    end
    return tempname()
  end
  return tempname(path; cleanup=cleanup)
end

# An Array where each element is written to disk
struct SerializedElementArray{T,N} <: AbstractArray{T,N}
  pathname::String
  dims::NTuple{N,Int}
  function SerializedElementArray{T,N}(
    ::UndefInitializer, d::NTuple{N,Integer}; cleanup=cleanup_default(), path=tempdir()
  ) where {T,N}
    mkpath(path)
    pathname = temppath(path; cleanup=cleanup)
    mkpath(pathname)
    return new{T,N}(pathname, d)
  end
end

const SerializedElementVector{T} = SerializedElementArray{T,1}
const SerializedElementMatrix{T} = SerializedElementArray{T,2}

pathname(d::SerializedElementArray) = d.pathname

function _filename(d::SerializedElementArray, I::Integer)
  return joinpath(d.pathname, "$I.bin")
end
filename(d::SerializedElementArray, I::Integer) = _filename(d, I)
filename(d::SerializedElementVector, I::Integer) = _filename(d, I)
function filename(d::SerializedElementArray, I::Integer...)
  return filename(d, CartesianIndex(I))
end
function filename(d::SerializedElementArray, I::CartesianIndex)
  return _filename(d, LinearIndices(d)[I])
end

Base.size(d::SerializedElementArray) = d.dims

# Check if the file exists, if not throw an UndefRefError
# Used in getindex and setindex! for SerializedElementArray
function checkfile(filename::AbstractString)
  !isfile(filename) && throw(UndefRefError())
  return nothing
end

function Base.getindex(d::SerializedElementArray{T}, I...)::T where {T}
  @boundscheck checkbounds(d, I...)
  filename_I = filename(d, I...)
  checkfile(filename_I)
  return deserialize(filename_I)
end

function Base.setindex!(d::SerializedElementArray, v, I...)
  @boundscheck checkbounds(d, I...)
  filename_I = filename(d, I...)
  serialize(filename_I, convert(eltype(d), v))
  return d
end

# Constructors with undefined values
function SerializedElementArray{T,N}(
  ::UndefInitializer, d::Vararg{Integer,N}; kw...
) where {T,N}
  return SerializedElementArray{T,N}(undef, d; kw...)
end
function SerializedElementArray{T,N}(; kw...) where {T,N}
  return SerializedElementArray{T,N}(undef, ntuple(_ -> 0, Val(N)); kw...)
end
function SerializedElementArray{T}(
  ::UndefInitializer, d::NTuple{N,Integer}; kw...
) where {T,N}
  return SerializedElementArray{T,N}(undef, d; kw...)
end
function SerializedElementArray{T}(
  ::UndefInitializer, d::Vararg{Integer,N}; kw...
) where {T,N}
  return SerializedElementArray{T,N}(undef, d; kw...)
end
function SerializedElementArray{<:Any,N}(
  ::UndefInitializer, d::NTuple{N,Integer}; kw...
) where {N}
  return SerializedElementArray{Any,N}(undef, d; kw...)
end
function SerializedElementArray{<:Any,N}(
  ::UndefInitializer, d::Vararg{Integer,N}; kw...
) where {N}
  return SerializedElementArray{Any,N}(undef, d; kw...)
end
function SerializedElementArray{<:Any,N}(; kw...) where {N}
  return SerializedElementArray{Any,N}(undef, ntuple(_ -> 0, Val(N)); kw...)
end

function SerializedElementArray{T,N}(A::AbstractArray; kw...) where {T,N}
  d = SerializedElementArray{T,N}(undef, size(A); kw...)
  for I in eachindex(A)
    if isassigned(A, I)
      d[I] = A[I]
    end
  end
  return d
end
function SerializedElementArray{T}(A::AbstractArray{<:Any,N}; kw...) where {T,N}
  return SerializedElementArray{T,N}(A; kw...)
end
function SerializedElementArray{<:Any,N}(A::AbstractArray{T}; kw...) where {T,N}
  return SerializedElementArray{T,N}(A; kw...)
end
function SerializedElementArray(A::AbstractArray{T,N}; kw...) where {T,N}
  return SerializedElementArray{T,N}(A; kw...)
end

"""
    SerializedElementArrays.disk(array::AbstractArray; cleanup=true, path=tempdir())

Convert the `AbstractArray` `array` to an array saved on disk of type `SerializedElementArrays.SerializedElementArray`. Each element of the array will be saved in serialized format to an individual file in a randomly generated directory.

If an `array` of type `SerializedElementArray` is input, this will simply return the original `array` (it acts like a conversion to type `SerializedElementArray`), ignoring the keyword arguments.

# Arguments
- `array::AbstractArray`: the array to convert to a `SerializedElementArrays.SerializedElementArray` stored on disk.

# Keywords
- `cleanup::Bool=true`: controls whether the process attempts to delete the returned path automatically when the process exits.
- `path=tempdir()`: the root of the path where a temporary directory will be made where the elements of the array will be saved. A directory with a random name assigned by the Base function `tempname` will be used.

!!! compat "Julia 1.4"
    This makes use of the Julia function `tempname` to create a temporary directory where the elements of the array will be saved. `tempname` only supports customizing the path and automatically cleaning up the files in Julia 1.4 and later. The `path` and `cleanup` arguments are only supported in Julia 1.4 and later, and in those previous versions they will default to `cleanup=false` and `path=tempdir()`.
"""
disk(array::AbstractArray; kw...) = SerializedElementArray(array; kw...)

disk(array::SerializedElementArray; kw...) = array

function disk(f::Function, args...; force_gc=true, full=false)
  d = disk(f(args...))
  force_gc && GC.gc(full)
  return d
end

end


text:
"""
    Sampler(y::Vector{Int}, X::Matrix{Float64}; kwargs...)

Initialize a Gibbs Sampler for the following Bayesian Negative-Binomial model:

```math
\\begin{aligned}
y_i | x_i, \\beta, \\gamma
&\\sim
\\text{NB}(r_{0y}, 1 / (1 + e^{- x_{i\\gamma}'\\beta_{\\gamma}})),
\\\\
\\beta_k | \\gamma_k
&\\sim
\\begin{cases}
    \\mathcal{N}(\\mu_{0\\beta}, \\Sigma_{0\\beta}),
    &\\text{if } \\gamma = 1,
    \\\\
    \\delta_0
    &\\text{if } \\gamma = 0,
\\end{cases}
\\\\
\\gamma
&\\sim
\\text{Womack}(\\eta_{0\\gamma}),
\\\\
s
&\\sim
\\text{Gamma}(a_{0s}, b_{0s}),
\\end{aligned}
```

given a response vector `y` and a design matrix `X`. This constructor does not
copy any of its arguments. Hence, for example, if `y[1]` is changed, it will
affect the Gibbs sampler.

where ``y_i \\in \\mathbb{N}_0`` is the response for the ``i``-observation,
``x_i \\in \\mathbb{R}^K`` is the associated covariate vector,
``\\delta_0(\\cdot)`` is the Dirac measure centred at ``0`` and
``\\text{Womack}(\\eta_0)`` is the Womack distribution with
parameter ``\\eta_{0\\gamma}`` on ``\\{0, 1\\}^K``.

# Keyword arguments

* `β = zeros(size(X, 2))`: current state of ``\\beta``.
* `Σ0β = 10 * I(size(X, 2))`: ``\\Sigma_{0\\beta}``.
* `μ0β = zeros(size(X, 2))`: ``\\mu_{0\\beta}``.
* `a0s = 1.0`: ``a_{0s}``.
* `b0s = 1.0`: ``a_{0s}``.

# Example

```julia
julia> using Random, Distributions
julia> rng = MersenneTwister(1)
julia> X = randn(rng, 100, 2)
julia> y = rand(rng, 0:2, 10)
julia> s = BayesNegativeBinomial.Sampler(y, X)
```
"""
struct Sampler
    y::Vector{Int}
    X::Matrix{Float64}
    mapping::Vector{Vector{Int}}
    β::Vector{Float64}
    ω::Vector{Float64}
    ξ::Vector{Float64}
    ϕ::Vector{Float64}
    q::Vector{Float64}
    g::Vector{Bool}
    A::Matrix{Float64}
    b::Vector{Float64}
    s::Vector{Int}
    a0s::Float64
    b0s::Float64
    μ0β::Vector{Float64}
    Σ0β::Matrix{Float64}
    ζ0g::Float64
    update_g::Vector{Bool}
    function Sampler(
        y::Vector{Int},
        X::Matrix{Float64};
        β::Vector{Float64} = zeros(size(X, 2)),
        s::Vector{Int} = [2],
        μ0β::Vector{Float64} = zeros(size(X, 2)),
        Σ0β::Matrix{Float64} = Matrix{Float64}(10 * I(size(X, 2))),
        mapping::Vector{Vector{Int}} = [[i] for i in 1:size(X, 2)],
        a0s::Float64 = 1.0,
        b0s::Float64 = 1.0,
        ζ0g::Float64 = 1.0,
        update_g = ones(Bool, length(mapping))
    )
        N, D = size(X)
        ω = zeros(N)
        ξ = zeros(N)
        ϕ = zeros(N)
        q = zeros(N)
        g = ones(Bool, length(mapping))
        A = zeros(D, D)
        b = zeros(D)
        new(y, X, mapping, β, ω, ξ, ϕ, q, g, A, b, s, a0s, b0s, μ0β, Σ0β, ζ0g, update_g)
    end
end

"""
    sample(rng::AbstractRNG, s::BayesNegativeBinomial.Sampler; kwargs...)

Draw a posterior sample using the Gibbs sampler `s`,
following Polson et al. (2013).

# Arguments

* `mcmcsize = 2000`: the total number of MCMC iterations.
* `burnin = 1000`: the number of *burn-in* iterations.

# Example

```julia
julia> using Random
julia> rng = MersenneTwister(1)
julia> X = randn(rng, 100, 2)
julia> y = rand(rng, 0:2, 10)
julia> sampler = BayesNegativeBinomial.Sampler(y, X)
julia> chain = BayesNegativeBinomial.sample(rng, sampler)
```

# References

1. Polson, N., Scott, J. & Windle, J. (2013) Bayesian inference for logistic
    models using Pólya–Gamma latent variables, *Journal of the American
    Statistical Association*, 108:504, 1339-1349,
    <https://doi.org/10.1080/01621459.2013.829001>.
"""
function sample(rng::AbstractRNG, sampler::Sampler; mcmcsize = 10000, burnin = 5000)
    chain = [zeros(size(sampler.X, 2)) for _ in 1:(mcmcsize - burnin)]
    for iter in 1:mcmcsize
        step!(rng, sampler)
        if iter > burnin
            chain[iter - burnin] .= sampler.β
        end
    end
    return chain
end

"""
    step!(rng::AbstractRNG, s::BayesNegativeBinomial.Sampler)

Perform 1 iteration of the Gibbs sampler `s`, following Polson et al. (2013).

# Example

```julia
julia> using Random, Distributions
julia> rng = MersenneTwister(1)
julia> X = randn(rng, 100, 2)
julia> y = rand(rng, 0:2, 10)
julia> s = BayesNegativeBinomial.Sampler(y, X)
julia> BayesNegativeBinomial.step!(rng, s)
```

# References

1. Polson, N., Scott, J. & Windle, J. (2013) Bayesian inference for logistic
    models using Pólya–Gamma latent variables, *Journal of the American
    Statistical Association*, 108:504, 1339-1349,
    <https://doi.org/10.1080/01621459.2013.829001>.
"""
function step!(rng::AbstractRNG, sampler::Sampler)
    step_ϕ!(sampler)
    step_ξ!(sampler)
    # step_s!(rng, sampler)
    step_ω!(rng, sampler)
    step_A!(sampler)
    step_b!(sampler)
    step_g!(rng, sampler)
    step_β!(rng, sampler)
    return nothing
end

function step_ξ!(sampler::Sampler)
    (; ξ, X, β) = sampler
    mul!(ξ, X, β)
    return nothing
end

function step_ϕ!(sampler::Sampler)
    (; ξ, ϕ) = sampler
    @. ϕ = 1.0 / (1.0 + exp(ξ))
    return nothing
end

function step_ω!(rng::AbstractRNG, sampler::Sampler)
    (; y, ω, ξ, s) = sampler
    for i in 1:length(ω)
        ω[i] = rand(rng, PolyaGammaPSWSampler(y[i] + s[], ξ[i]))
    end
    return nothing
end

function step_g!(rng::AbstractRNG, sampler::Sampler)
    (; update_g, mapping, g, μ0β, Σ0β, ζ0g) = sampler
    D = length(μ0β)
    gexp = zeros(Bool, D)
    for d in 1:length(mapping)
        gexp[mapping[d]] .= g[d]
    end
    pg = Womack(length(g), ζ0g)
    for d in 1:length(g)
        update_g[d] || continue
        logodds = 0.0
        for val in 0:1
            g[d] = val
            gexp[mapping[d]] .= val
            m1, Σ1 = posterior_hyperparameters(sampler)
            logodds += (-1)^(val + 1) * (
                logpdf(pg, g) +
                logpdf(MvNormal(μ0β[gexp], Σ0β[gexp, gexp]), zeros(sum(gexp))) -
                logpdf(MvNormal(m1, Σ1), zeros(length(m1)))
            )
        end
        g[d] = rand(rng) < exp(logodds) / (1.0 + exp(logodds))
        gexp[mapping[d]] .= g[d]
    end
    return nothing
end

function step_β!(rng::AbstractRNG, sampler::Sampler)
    (; mapping, β, g) = sampler
    D = length(β)
    gexp = zeros(Bool, D)
    for d in 1:length(mapping)
        gexp[mapping[d]] .= g[d]
    end
    β .= 0.0
    m1, Σ1 = posterior_hyperparameters(sampler)
    β[gexp] .= rand(rng, MvNormal(m1, Σ1))
    return nothing
end

function posterior_hyperparameters(sampler::Sampler)
    (; mapping, g, A, b, μ0β, Σ0β) = sampler
    D = length(μ0β)
    gexp = zeros(Bool, D)
    for d in 1:length(mapping)
        gexp[mapping[d]] .= g[d]
    end
    Σ1 = inv(cholesky(Symmetric(A[gexp, gexp])))
    m1 = Σ1 * (b[gexp] + Σ0β[gexp, gexp] \ μ0β[gexp])
    return m1, Σ1
end

function step_A!(sampler::Sampler)
    (; X, ω, Σ0β, A) = sampler
    A .= X' * Diagonal(ω) * X + inv(Σ0β)
    return nothing
end

function step_b!(sampler::Sampler)
    (; X, y, s, b) = sampler
    b .= X' * (y .- s[]) / 2
    return nothing
end

# function step_s!(rng, sampler::Sampler)
#     (; y, q, ϕ, s, a0s, b0s) = sampler
#     N = length(q)
#     s0 = s[]
#     for i in 1:N
#         q[i] = 0
#         for j in 1:y[i]
#             q[i] += rand(rng) <= s0 / (s0 + j - 1)
#         end
#     end
#     a1s = a0s + sum(q)
#     b1s = b0s
#     for i in 1:N
#         b1s - log(1 - ϕ[i])
#     end
#     ds = Gamma(a1s, 1 / b1s)
#     s[] = rand(rng, ds)
#     return nothing
# end

struct Womack <: DiscreteMultivariateDistribution
    D::Int
    ζ::Float64
    p::Vector{Float64}
    punnormalized::Vector{Float64}
    function Womack(D::Int, ζ::Float64)
        p = big.([zeros(D); 1.0])
        for d1 in (D - 1):-1:0
            for d2 in 1:(D - d1)
                p[1 + d1] += ζ * p[1 + d1 + d2] * binomial(big(d1 + d2), big(d1))
            end
        end
        p /= sum(p)
        punnormalized = copy(p)
        for d1 in 1:D
            p[d1] /= binomial(big(D), big(d1 - 1))
        end
        return new(D, ζ, p, punnormalized)
    end
end

function pdf(d::Womack, g::Vector{Bool})
    return d.p[sum(g) + 1]
end

function logpdf(d::Womack, g::Vector{Bool})
    return log(pdf(d, g))
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
using Base.BinaryPlatforms

# copied from libsingular_julia:
# See https://github.com/JuliaLang/Pkg.jl/issues/2942
# Once this Pkg issue is resolved, this must be removed
uuid = Base.UUID("a83860b7-747b-57cf-bf1f-3e79990d037f")
delete!(Pkg.Types.get_last_stdlibs(v"1.6.3"), uuid)

name = "libpolymake_julia"
version = v"0.8.0"

julia_versions = [v"1.6.3", v"1.7.0", v"1.8.0"]

# Collection of sources required to build libpolymake_julia
sources = [
    ArchiveSource("https://github.com/oscar-system/libpolymake-julia/archive/v$(version).tar.gz",
                  "da459c1fc819a446cfa683cc8e6e890a54194b7a943ac8d62c4c71061512dbee"),
]

# Bash recipe for building across all platforms
script = raw"""
# change default perl which interferes with the hostbuild perl
rm -f /usr/bin/perl
ln -s $host_bindir/perl /usr/bin/perl

cmake libpolymake-j*/ -B build \
   -DJulia_PREFIX="$prefix" \
   -DCMAKE_INSTALL_PREFIX="$prefix" \
   -DCMAKE_FIND_ROOT_PATH="$prefix" \
   -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
   -DCMAKE_BUILD_TYPE=Release

VERBOSE=ON cmake --build build --config Release --target install -- -j${nproc}

install_license libpolymake-j*/LICENSE.md

jsondir=${prefix}/share/libpolymake_julia/appsjson/
mkdir -p $jsondir
$host_bindir/perl $host_bindir/polymake --iscript libpolymake-j*/src/polymake/apptojson.pl $jsondir
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")

platforms = vcat(libjulia_platforms.(julia_versions)...)
filter!(p -> !Sys.iswindows(p) && arch(p) != "armv6l", platforms)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libpolymake_julia", :libpolymake_julia),
    FileProduct("share/libpolymake_julia/type_translator.jl",:type_translator),
    FileProduct("share/libpolymake_julia/appsjson",:appsjson),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("libjulia_jll"),
    BuildDependency("GMP_jll"),
    BuildDependency("MPFR_jll"),
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("FLINT_jll", compat = "~200.800.401"),
    Dependency("TOPCOM_jll"),
    Dependency("lib4ti2_jll"),
    Dependency("libcxxwrap_julia_jll"),
    Dependency("polymake_jll"; compat = "~400.600.0"),

    HostBuildDependency(PackageSpec(name="Perl_jll", version=v"5.34.0")),
    HostBuildDependency(PackageSpec(name="polymake_jll", version=v"400.600.0")),
    HostBuildDependency("lib4ti2_jll"),
    HostBuildDependency("TOPCOM_jll"),
]

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

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

if !((@__DIR__) ∈ LOAD_PATH)
    push!(LOAD_PATH, (@__DIR__))
end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
function Base.convert(T::Type{Vector{SVector{N2,T1}}}, M::Matrix{T2}) where {N2,T1,T2}
    N1 = size(M)[1]
    @assert size(M)[2] == N2
    Mout = [@SVector zeros(T1, N2) for i = 1:N1]
    for i = 1:N1
        Mout[i] = convert(SVector{N2,T1}, M[i,:])
    end
    return Mout
end

function deleteat(M::Array{T,2}, i1::Int64, i2::Int64) where T
    [M[1:i1 - 1,1:i2 - 1] M[1:i1 - 1,i2 + 1:end];M[i1 + 1:end,1:i2 - 1] M[i1 + 1:end,i2 + 1:end]]
end

deleteat(M::Array{T,2},i::Int64) where T = deleteat(M, i, i)

skew(v::AbstractVector{T}) where T = SMatrix{3,3,T,9}(0, v[3], -v[2], -v[3], 0, v[1], v[2], -v[1], 0)

Base.:*(u::LinearAlgebra.AdjointAbsVec, v::SVector{1,T}) where T = u * v[1]
Base.:*(u::AbstractVector, v::SVector{1,T}) where T = u * v[1]

@inline svcat(a::T) where T = SVector{1,T}(a)
@inline svcat(a::StaticArray) = a
@inline svcat(a::T, b::T) where T = SVector{2,T}(a, b)
@inline svcat(a::StaticArray, b::StaticArray) = vcat(a, b)
@inline svcat(a::StaticArray{Tuple{N},T,1}, b::T) where {T,N} = vcat(a, SVector{1,T}(b))
@inline svcat(a::T, b::StaticArray{Tuple{N},T,1}) where {T,N} = vcat(SVector{1,T}(a), b)

@inline svcat(a, b, c...) = svcat(svcat(a, b), svcat(c...))


text:
using ELM

elm = ExtremeLearningMachine(100)

# Testing XOR
x = [1.0 1.0; 0.0 1.0; 0.0 0.0; 1.0 0.0]
y = [0.0, 1.0, 0.0, 1.0]

fit!(elm, x, y)

y_pred = predict(elm, [1 1.; 0 1; 1 1])
@assert y_pred[1] < 0.2
@assert y_pred[2] > 0.8
@assert y_pred[3] < 0.2


text:
## Test π₀ estimator methods ##
module Test_pi0

using MultipleTesting
using Test
using StatsBase

include("utils.jl")

@testset "π₀ estimators" begin

    ## test case: deterministic
    p0 = collect(0.01:0.01:1)
    p1 = p0.^10
    p = [p0; p1] ## unsorted
    pi0 = length(p0) / length(p)

    lambdas = collect(0.05:0.05:0.95)

    @testset "Storey π₀" begin

        @test typeof(Storey()) <: Pi0Estimator
        @test typeof(Storey(0.5)) <: Pi0Estimator
        @test estimate(p, Storey(0.2)) ≈ 0.6
        @test estimate(p, Storey(0.0)) ≈ 1.0
        @test estimate(p, Storey(1.0)) ≈ 1.0

        p_unsort = unsort(p)
        @test !issorted(p_unsort)
        @test estimate(p_unsort, Storey(0.2)) ≈ 0.6
        @test !issorted(p_unsort)

        @test_throws DomainError Storey(-0.1)
        @test_throws DomainError Storey(1.1)

    end

    @testset "StoreyBootstrap π₀" begin

        q = 0.1

        @test typeof(StoreyBootstrap()) <: Pi0Estimator

        @test estimate(p, StoreyBootstrap()) ≈ 0.6
        @test estimate(p0, StoreyBootstrap()) ≈ 1.0
        @test estimate(p1, StoreyBootstrap()) ≈ 0.15

        @test estimate(p, StoreyBootstrap(lambdas, q)) ≈ 0.6
        @test estimate(p0, StoreyBootstrap(lambdas, q)) ≈ 1.0
        @test estimate(p1, StoreyBootstrap(lambdas, q)) ≈ 0.15

        ## unsorted lambdas
        lambdas_unsort = unsort(lambdas)
        @test !issorted(lambdas_unsort)
        @test estimate(p, StoreyBootstrap(lambdas_unsort, q)) ≈ 0.6
        @test !issorted(lambdas_unsort)

        @test_throws DomainError StoreyBootstrap(lambdas, -0.1)
        @test_throws DomainError StoreyBootstrap(lambdas, 1.1)
        @test_throws MethodError StoreyBootstrap(0.5)
        @test_throws MethodError StoreyBootstrap(lambdas)
        @test_throws MethodError StoreyBootstrap(0.5, lambdas)

    end

    @testset "LeastSlope π₀" begin

        # alternative, vectorized version
        # used for comparison and compactness
        function lsl_pi0_vec(pValues::AbstractVector{T}) where T <: AbstractFloat
            n = length(pValues)
            pValues = MultipleTesting.sort(pValues)
            s = (1 .- pValues) ./ (n:-1:1)
            d = diff(s) .< 0
            idx = something(findfirst(d), 0) + 1
            pi0 = min(1 / s[idx] + 1, n) / n
            return pi0
        end

        @test typeof(LeastSlope()) <: Pi0Estimator

        ## checked against structSSI::pi0.lsl
        @test isapprox(estimate(p, LeastSlope()), 0.62, atol = 1e-2)
        @test isapprox(estimate(p0, LeastSlope()), 1.0, atol = 1e-2)
        @test isapprox(estimate(p1, LeastSlope()), 0.16, atol = 1e-2)

        ## check against internal reference implementation
        @test estimate(p, LeastSlope()) ≈ lsl_pi0_vec(p)
        @test estimate(p0, LeastSlope()) ≈ lsl_pi0_vec(p0)
        @test estimate(p1, LeastSlope()) ≈ lsl_pi0_vec(p1)

        @test_throws MethodError LeastSlope(0.1)

        ## unsorted p-values
        p_unsort = unsort(p)
        @test !issorted(p_unsort)
        @test isapprox(estimate(p_unsort, LeastSlope()), 0.62, atol = 1e-2)
        @test !issorted(p_unsort)

        p_unsort = unsort(p)
        @test !issorted(p_unsort)
        @test isapprox(lsl_pi0_vec(p_unsort), 0.62, atol = 1e-2)
        @test !issorted(p_unsort)

    end

    @testset "Oracle π₀" begin

        @test estimate(p, Oracle(0.5)) == 0.5
        @test estimate(p0, Oracle(0.6)) == 0.6
        @test estimate(p1, Oracle()) == 1.0

    end

    @testset "TwoStep π₀" begin

        alpha = 0.05

        @test typeof(TwoStep(alpha)) <: Pi0Estimator

        ## checked against mutoss::TSBKY_pi0_est
        @test estimate(p, TwoStep()) ≈ 0.665
        @test estimate(p, TwoStep(alpha)) ≈ 0.665
        @test estimate(p0, TwoStep(alpha)) ≈ 1.0
        @test estimate(p1, TwoStep(alpha)) ≈ 0.29
        @test estimate(p, TwoStep(alpha, BenjaminiHochberg())) ≈ 0.665

        @test estimate(p, TwoStep(0.1)) ≈ 0.63
        @test estimate(p, TwoStep(0.0)) ≈ 1.0
        @test estimate(p, TwoStep(1.0)) ≈ 0.415
        @test estimate(p, TwoStep(0.1, BenjaminiHochberg())) ≈ 0.63

        ## unsorted p-values
        p_unsort = unsort(p)
        @test !issorted(p_unsort)
        @test estimate(p_unsort, TwoStep(alpha)) ≈ 0.665
        @test !issorted(p_unsort)

    end

    @testset "RightBoundary π₀" begin

        # R code tested against (note this uses equidistant grid so we have to use
        # λseq as in Storey for comparison
        # ```R
        # library(pi0)
        # p0 <- seq(0.01,1, by=0.01)
        # histf1(p0, rightBoundary=TRUE)
        # [1] 1
        # p1 <- p0^10
        # histf1(p1, rightBoundary=TRUE)
        # [1] 0.1428571
        # p <- c(p0,p1)# histf1(p, rightBoundary=TRUE)
        # [1] 0.5714286
        # ```

        @test typeof(RightBoundary()) <: Pi0Estimator
        @test typeof(RightBoundary(lambdas)) <: Pi0Estimator

        # only eps because pi0 package uses right closed histograms
        @test isapprox(estimate(p, RightBoundary(lambdas)), 0.5714286, atol = 0.02)
        @test estimate(p0, RightBoundary(lambdas)) ≈ 1.0
        @test isapprox(estimate(p1, RightBoundary(lambdas)), 0.1428571, atol = 1e-7)

        @test estimate(p, RightBoundary(lambdas)) ≈ estimate(p, RightBoundary(unsort(lambdas)))

        # not checked against R implementation but should hold (check for default lambda grid)
        @test estimate(p0, RightBoundary()) ≈ 1.0

        p_unsort = unsort(p1)
        @test !issorted(p_unsort)
        @test isapprox(estimate(p_unsort, RightBoundary(lambdas)), 0.1428571, atol = 1e-7)
        @test !issorted(p_unsort)

    end

    @testset "CensoredBUM π₀" begin

        function cbum_pi0_naive(pValues::AbstractVector{T},
                γ0::AbstractFloat = 0.5, λ::AbstractFloat = 0.05,
                xtol::AbstractFloat = 1e-6, maxiter::Integer = 10000) where T <: AbstractFloat

            n = length(pValues)
            z = fill(1 - γ0, n)
            idx_left = pValues .< λ
            idx_right = .!idx_left
            pi0_old = γ0 = α = γ = Inf
            lpr = log.(pValues[idx_right])
            ll = log(λ)
            for i in 1:maxiter
            γ = sum(1 .- z) / n
            α = -sum(z[idx_right])
            α = α / ( ll * sum(z[idx_left]) + sum(z[idx_right] .* lpr) )
            xl = (1 - γ) * (λ^α)
            z[idx_left] .= xl / (γ * λ + xl)
            xr = (1 - γ) * α * pValues[idx_right].^(α - 1)
            z[idx_right] = xr ./ (γ .+ xr)
            pi0_new = γ + (1 - γ) * α
            if abs(pi0_new - pi0_old) <= xtol
            return pi0_new, [γ, α], true
            end
            γ0 = γ
            pi0_old = pi0_new
            end
            return NaN, [γ, α], false
        end

        @test typeof(CensoredBUM()) <: Pi0Estimator
        @test typeof(CensoredBUM(0.2, 0.1)) <: Pi0Estimator

        @test isapprox(estimate(p, CensoredBUM()), 0.55797, atol = 1e-5)
        @test isapprox(estimate(p0, CensoredBUM()), 1.0, atol = 1e-5)
        @test isapprox(estimate(p1, CensoredBUM()), 0.11608, atol = 2e-5)

        # test against internal reference implementation
        # p0 case is not handled well by naive implementation
        @test estimate(p, CensoredBUM()) ≈ cbum_pi0_naive(p)[1]
        @test estimate(p1, CensoredBUM()) ≈ cbum_pi0_naive(p1)[1]

        ## test case that does not converge
        @test isnan(estimate(p, CensoredBUM(0.5, 0.05, 1e-6, 2)))

        @test_throws DomainError CensoredBUM(-0.5, 0.05)
        @test_throws DomainError CensoredBUM(1.5, 0.05)
        @test_throws DomainError CensoredBUM(0.5, -0.05)
        @test_throws DomainError CensoredBUM(0.5, 1.05)

        @test_throws DomainError CensoredBUM(0.5, 0.05, -1e-6, 100)
        @test_throws DomainError CensoredBUM(0.5, 0.05, 1.1, 100)
        @test_throws DomainError CensoredBUM(0.5, 0.05, 1e-6, -10)

        f = fit(CensoredBUM(), p)
        @test typeof(f) <: CensoredBUMFit
        @test isapprox(f.π0, 0.55797, atol = 1e-5)

        # denominator becomes 0 if all p-values are 1
        @test isnan(estimate(ones(20), CensoredBUM()))

        pi0_est, pars, is_converged = MultipleTesting.cbum_pi0(ones(20))
        @test isnan(pi0_est)
        @test !is_converged

        pi0_est, pars, is_converged = cbum_pi0_naive(ones(20))
        @test isnan(pi0_est)
        @test !is_converged

    end

    @testset "BUM π₀" begin

        @test isapprox(estimate(p, BUM(0.5)), 0.55528, atol = 1e-5)
        @test isapprox(estimate(p, BUM()), 0.55528, atol = 1e-5)
        @test isapprox(estimate(p0, BUM()), 1.0, atol = 1e-5)
        @test isapprox(estimate(p1, BUM()), 0.10874, atol = 1e-5)

        ## test case that does not converge
        @test isnan(estimate(p, BUM(0.5, 1e-6, 2)))

        @test typeof(BUM()) <: Pi0Estimator

        @test_throws DomainError BUM(-0.5)
        @test_throws DomainError BUM(1.5)

    end

    @testset "FlatGrenander π₀" begin

        @test typeof(FlatGrenander()) <: Pi0Estimator

        pu = collect(0.1:0.05:0.9)
        @test estimate(pu, FlatGrenander()) ≈ 1.0
        @test estimate(pu.^0.5, FlatGrenander()) ≈ 1.0

        # no reference values from literature/software available
        @test estimate(p0, FlatGrenander()) ≈ 1.0
        @test isapprox(estimate(p1, FlatGrenander()), 0.10458, atol = 1e-4)
        @test isapprox(estimate(p, FlatGrenander()), pi0, atol = 0.1)

        ## longest constant interval: low level
        lci = MultipleTesting.longest_constant_interval

        px = [0:0.1:1;];
        f = [1.5, 1.5, 1.5, 1.5, 1.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.1];
        @test lci(px, f) ≈ 0.5

        f = [1.5, 1.5, 1.5, 1.5, 1.5, 1.2, 1.2, 1.2, 0.2, 0.2, 0.1];
        @test lci(px, f) ≈ 0.2

        f = [0.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.1];
        @test lci(px, f) ≈ 0.2

        f = [0.5, 0.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.2, 0.2, 0.2, 0.1];
        @test lci(px, f) ≈ 0.5

        px = [0.1, 0.3, 0.5, 0.9];
        f = [0.5, 0.5, 0.2, 0.2];
        @test lci(px, f) ≈ 0.2

        px = [0.1, 0.5, 0.7, 0.9];
        f = [0.5, 0.5, 0.2, 0.2];
        @test lci(px, f) ≈ 0.5

    end

    @testset "ConvexDecreasing π₀" begin

        @test typeof(ConvexDecreasing()) <: Pi0Estimator
        @test typeof(ConvexDecreasing(100, 1e-6, 1000)) <: Pi0Estimator

        @test isapprox(estimate(p, ConvexDecreasing()), 0.5739054, atol = 1e-4)
        @test isapprox(estimate(p0, ConvexDecreasing()), 1.0, atol = 1e-6)
        @test isapprox(estimate(p1, ConvexDecreasing()), 0.1323284, atol = 1e-3)

        @test isapprox(estimate(ones(20), ConvexDecreasing()), 1.0, atol = 1e-6)

        ## test case that does not converge
        @test isnan(estimate(p, ConvexDecreasing(100, 1e-6, 2)))

        @test_throws DomainError ConvexDecreasing(100, 2.0, 1000)
        @test_throws DomainError ConvexDecreasing(100, -1.0, 1000)
        @test_throws DomainError ConvexDecreasing(-100, 0.01, 1000)
        @test_throws DomainError ConvexDecreasing(100, 0.01, -1000)

        f = fit(ConvexDecreasing(), p)
        @test typeof(f) <: ConvexDecreasingFit
        @test isapprox(f.π0, 0.5739054, atol = 1e-4)

    end

end

end


text:
using IntervalArithmetic
if VERSION < v"0.7.0-DEV.2004"
    using Base.Test
else
    using Test
end

@testset "Complex interval operations" begin
    a = @interval 1im
    @test typeof(a)== Complex{IntervalArithmetic.Interval{Float64}}
    @test a ==  Interval(0) + Interval(1)*im
    @test a * a == Interval(-1)
    @test a + a == Interval(2)*im
    @test a - a == 0
    @test a / a == 1
    @test a^2 == -1
end

@testset "Complex functions" begin
    Z = (3 ± 1e-7) + (4 ± 1e-7)*im
    @test sin(Z) == complex(sin(real(Z))*cosh(imag(Z)),sinh(imag(Z))*cos(real(Z)))
    @test exp(-im * Interval(π)) == Interval(-1.0, -0.9999999999999999) - Interval(1.224646799147353e-16, 1.2246467991473532e-16)*im
end


text:
using compost
using SafeTestsets
using DataFrames
using Test

# Test functions that calculate comparative advantage: ==============
@safetestset "Comparative advantage tests" begin include(normpath(@__DIR__, "test/comparative-advantage-tests.jl")) end

# Test functions the calculate similarity and distance measures: ====
@safetestset "Relatedness tests" begin include(normpath(@__DIR__, "test/relatedness-tests.jl")) end

# Test complexity algorithms: =======================================
@safetestset "Complexity algorithm tests" begin include(normpath(@__DIR__, "test/similarity-distance-tests.jl")) end

# Test forecasting functions: =======================================
# TODO:

max_stars_repo_name: sbpost/compost

max_issues_repo_name: sbpost/compost

max_forks_repo_name: sbpost/compost


text:
using Convex, SCS

function EMAC(Y,N,K,omega,sigma)
# Nuclear norm minimization for line spectral estimation from measurements
# of a sparse linear array
#
# Inputs:
# Y: Data matrix
# N: ambient dimension \geq M
# K: number of sinusoids
# omega: indices of the observed entries of the data in the ambient space
# sigma: noise variance of the additive white Gaussian noise with mean zero
#
# M = size(Y,1): number of measurements per snapshot
# L = size(Y,2): number of snapshots
#
# Output:
# f_hat: estimate of the frequencies of the K sinusoids in [rad]

(M,L) = size(Y)
m::Int = ceil(N/2)
n = N + 1 - m
lambda = sigma*L*length(omega)/(1e2*sqrt(N))
#lambda = 1e-4

Q1 = ComplexVariable(n*L,n*L)
Q2 = ComplexVariable(m,m)
H = ComplexVariable(m,n*L)
z = ComplexVariable(N,L)

constraints = ([Q1 H'; H Q2] in :SDP)
constraints += [H[i,j+k*n] == z[i+j-1,k+1] for i = 1:m, j = 1:n, k=0:L-1]
constraints += [sumsquares(abs(vec(z[omega,:]-Y))) < lambda]
objective = real( tr(Q1) + tr(Q2) )
problem = minimize(objective, constraints)

solve!(problem, SCS.Optimizer)

clmnsY0 = setdiff(1:n*L,(1:L).*n)
clmnsY1 = clmnsY0.+1

Y0 = H.value[:,clmnsY0]
Y1 = H.value[:,clmnsY1]

F = svd(Y0)
p = sortperm(F.S,rev=true)[1:K]
f = eigvals(Diagonal(ones(K)./F.S[p])*F.U[:,p]'*Y1*F.V[:,p])
f_hat = (angle.(f)/(2*pi) .+ 1).%1

return f_hat
end