text:
######################
####Utility Method####
######################
function attachSolverWithAttributesToJuMPModel(model:: Model, options::Dict{String, Any})
	set_optimizer(model, OnePhase.OnePhaseSolver)
	for (name, value) in options
        sname = string(name)
	set_optimizer_attribute(model, sname, value)
    end
end

######################
##### ROSENBROOK #####
######################
function rosenbrook1()
    model = Model()
    @variable(model, x)
    @variable(model, y)
    @NLobjective(model, Min, (2.0 - x)^2 + 100 * (y - x^2)^2)
    return model
end

#function test_rosenbrook1(solver)
function test_rosenbrook1(options::Dict{String, Any})
    @testset "test_rosenbrook1" begin
        model = rosenbrook1()
		attachSolverWithAttributesToJuMPModel(model, options)
        @test_broken optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test_broken status == :Optimal
    end
end

function rosenbrook2()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @NLobjective(model, Min, (2.0 - x)^2 + 100 * (y - x^2)^2)
    @constraint(model, x + y >= 0.1)
	@NLconstraint(model, x* y + x >= 0.1)
    return model
end

#function test_rosenbrook2(solver)
function test_rosenbrook2(options::Dict{String, Any})
    @testset "test_rosenbrook2" begin
        model = rosenbrook2()
		attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test status == :Optimal
        check_rosenbrook(model)
    end
end

function rosenbrook3()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @NLobjective(model, Min, (2.0 - x)^2 + 100 * (y - x^2)^2)
    @NLconstraint(model, x^2 + y^2 >= 0.5)
    return model
end

#function test_rosenbrook3(solver)
function test_rosenbrook3(options::Dict{String, Any})
    @testset "test_rosenbrook3" begin
        model = rosenbrook3()
		attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test status == :Optimal
        check_rosenbrook(model)
    end
end

function rosenbrook4()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    #FIXME
    @NLconstraint(model, (x + y) ^ 2 >= 0)
    @NLobjective(model, Min, (2.0 - x)^2 + 100 * (y - x^2)^2)
    return model
end

#function test_rosenbrook4(solver)
function test_rosenbrook4(options::Dict{String, Any})
    @testset "test_rosenbrook4" begin
        model = rosenbrook4()
		attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test status == :Optimal
    end
end

function check_rosenbrook(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) - 2.0) < tol
    @test abs(JuMP.value.(model[:y]) - 4.0) < tol
end

########################
##### FEASIBLE LPS #####
########################
function toy_lp0()
    model = Model()
    @variable(model, x)
    @objective(model, Min, x)
    @NLconstraint(model, x >= 4.0)
    return model
end

function toy_lp1()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @NLobjective(model, Min, -x - 100 * y)
    @constraint(model, x + y <= 1.0)
    return model
end

function check_toy_lp1(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) - 0.0) < tol
    @test abs(JuMP.value.(model[:y]) - 1.0) < tol
end

function test_toy_lp1(options::Dict{String, Any})
    model = toy_lp1()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    check_toy_lp1(model)
end

function toy_lp2()
    model = Model()
    @variable(model, 0.0 <= x <= 1.0)
    @variable(model, 0.0 <= y <= 1.0)
    @NLobjective(model, Min, -x - 100 * y)
    @constraint(model, x + y <= 2.0)
    return model
end

function check_toy_lp2(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) - 1.0) < tol
    @test abs(JuMP.value.(model[:y]) - 1.0) < tol
end

#function test_toy_lp2(solver)
function test_toy_lp2(options::Dict{String, Any})
    model = toy_lp2()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    check_toy_lp2(model)
end

function toy_lp3()
    model = Model()
    @variable(model, 0.0 <= x <= 1.0)
    @variable(model, 0.0 <= y <= 1.0)
    @NLobjective(model, Min, x)
    @constraint(model, 1.0 <= x + y <= 2.0)
    return model
end

function check_toy_lp3(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) - 0.0) < tol
    @test abs(JuMP.value.(model[:y]) - 1.0) < tol
end

#function test_toy_lp3(solver)
function test_toy_lp3(options::Dict{String, Any})
    model = toy_lp3()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    check_toy_lp3(model)
end

function check_toy_lp4(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) - 0.0) < tol
    @test abs(JuMP.value.(model[:y]) - 1.0) < tol
end

function toy_lp4()
    model = Model()
    @variable(model, x, lower_bound=0.0, upper_bound=1.0)
    @variable(model, y, lower_bound=0.0, upper_bound=1.0)
    @NLobjective(model, Min, x)
    @constraint(model, 1.0 <= x + y <= 2.0)
    return model
end

function check_toy_lp4(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) - 0.0) < tol
    @test abs(JuMP.value.(model[:y]) - 1.0) < tol
end

#function test_toy_lp4(solver)
function test_toy_lp4(options::Dict{String, Any})
    model = toy_lp4()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    check_toy_lp4(model)
end
#Order =, >, <, < <
function toy_lp5()
    model = Model()
	@variable(model, x, lower_bound=0.0, upper_bound=1.0)
    @variable(model, y, lower_bound=0.0, upper_bound=1.0)
    @NLobjective(model, Min, x)
    @constraint(model, x + y == 1.0)
    @constraint(model, x * 32.5 + y * 32.5 == 32.5)
    @constraint(model, 3.0 * x + 3.0 * y <= 3.0)
    return model
end

function test_toy_lp5(options::Dict{String, Any})
    model = toy_lp5()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    check_toy_lp4(model)
end

function toy_lp6()
    model = Model()
    @variable(model, x, lower_bound=0.0, upper_bound=1.0)
    @variable(model, y, lower_bound=0.0, upper_bound=1.0)
    @NLobjective(model, Min, x)
    @constraint(model, x + y == 1.0)
    @constraint(model, x * 5.5 + y * 5.5 == 5.5)
    return model
end

function test_toy_lp6(options::Dict{String, Any})
    model = toy_lp6()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    check_toy_lp4(model)
end

function toy_lp7()
    model = Model()
    @variable(model, x, lower_bound=0.0, upper_bound=1.0)
    @variable(model, y, lower_bound=0.0, upper_bound=1.0)
    @NLobjective(model, Min, x)
    @constraint(model, 2.0 * x + y == 1.0)
    return model
end

function test_toy_lp7(options::Dict{String, Any})
    model = toy_lp7()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    check_toy_lp4(model)
end

function toy_lp8()
    model = Model()
    @variable(model, x, lower_bound=0.0, upper_bound=1.0)
    @variable(model, y, lower_bound=0.0, upper_bound=1.0)
    @NLobjective(model, Min, x)
    @constraint(model, x + y >= 1.0)
    @constraint(model, x * 5.5 + y * 5.5 <= 5.5)
    return model
end

function test_toy_lp8(options::Dict{String, Any})
    model = toy_lp8()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    check_toy_lp4(model)
end

##########################
##### INFEASIBLE LPS #####
##########################

function toy_lp_inf1()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @NLobjective(model, Min, x + 100 * y)
    @constraint(model, x + 2 * y <= -1.0)
    return model
end

function toy_lp_inf2()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @NLobjective(model, Min, x + 100 * y)
    @constraint(model, x + 2 * y <= 2.0)
    @constraint(model, x + 2 * y >= 4.0)
    return model
end

#######################
##### CONVEX NLPS #####
#######################

function circle1()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @NLobjective(model, Min, x + 100 * y)
    @NLconstraint(model, x^2 + y^2 <= 1.0)
    @NLconstraint(model, (x-2.0)^2 + y^2 <= 1.0)
    return model
end

function check_circle1(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) - 1.0) < tol
    @test abs(JuMP.value.(model[:y]) - 0.0) < tol
end

function circle2()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @NLobjective(model, Min, x^3 + y^3)
    @NLconstraint(model, x^2 + y^2 <= 1.0)
    return model
end

function check_circle2(model)
    tol = 1e-2
    @test abs(JuMP.value.(model[:x]) - 0.0) < tol
    @test abs(JuMP.value.(model[:y]) - 0.0) < tol
end

function quad_opt()
    model = Model()
    @variable(model, x)
    @variable(model, y)
    @NLobjective(model, Min, y)
    @NLconstraint(model, x^2 <= y)
    return model
end

function check_quad_opt(model)
    tol = 1e-2
    @test abs(JuMP.value.(model[:x]) - 0.0) < tol
    @test abs(JuMP.value.(model[:y]) - 0.0) < tol
end

##########################
##### NONCONVEX NLPS #####
##########################
function circle_nc1()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @NLobjective(model, Min, x + 100 * y)
    @NLconstraint(model, x^2 + y^2 == 1.0)
    @NLconstraint(model, (x-2.0)^2 + y^2 == 1.0)
    return model
end

function check_circle_nc1(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) - 1.0) < tol
    @test abs(JuMP.value.(model[:y]) - 0.0) < tol
end

function circle_nc2()
    model = Model()
    @variable(model, x, start = 1.0)
    @variable(model, y, start = 1.0)
    @NLobjective(model, Min, x)
    @NLconstraint(model, x^2 + y^2 == 1.0)
    return model
end

function check_circle_nc2(model)
    tol = 1e-3
    @test abs(JuMP.value.(model[:x]) + 1.0) < tol
    @test abs(JuMP.value.(model[:y]) - 0.0) < tol
end

function circle_nc_inf1()
    model = Model()
    @variable(model, x, start = 1.0)
    @variable(model, y, start = 1.0)
    @NLobjective(model, Min, x)
    @NLconstraint(model, x^2 + y^2 == 1.0)
    @NLconstraint(model, x^2 + 2 * y^2 == 4.0)
    return model
end

##############################
##### UNBOUNDED PROBLEMS #####
##############################

function lp_unbd()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y)
    @NLobjective(model, Min, -x)
    @NLconstraint(model, x - y <= 1.0)
    return model
end

function circle_nc_unbd()
    model = Model()
    @variable(model, x)
    @variable(model, y)
    @NLobjective(model, Min, x + 0.1 * y)
    @NLconstraint(model, x^2 + y^2 >= 1.0)
    return model
end

function quad_unbd()
    model = Model()
    @variable(model, x)
    @variable(model, y)
    @NLobjective(model, Min, x)
    @NLconstraint(model, x^2 <= y)
    return model
end

#####################################
##### UNBOUNDED FEASIBLE REGION #####
#####################################

function unbd_feas()
    model = Model()
    @variable(model, x >= 0.0)
    @variable(model, y >= 0.0)
    @variable(model, z >= 0.0)
    @NLobjective(model, Min, y)
    @NLconstraint(model, x^2 <= y)
    @NLconstraint(model, z >= 0.0)
    return model
end

function test_unbd_feas(options::Dict{String, Any})
    println("test_unbd_feas")
    model = unbd_feas()
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal
    @test JuMP.value.(model[:z]) < 1e5
end

###########################
##### STARTING POINTS #####
###########################

function starting_point_prob(start::Float64)
    model = Model()
    @variable(model, x, start = start)
    @NLobjective(model, Min, -x^2)
    @NLconstraint(model, -1.0 <= x <= 1.0)
    return model
end

function test_starting_point(options::Dict{String, Any},starting_point::Float64)
    if starting_point == 0.0
        warn("don't select this as a starting point")
    end
    model = starting_point_prob(starting_point)
    attachSolverWithAttributesToJuMPModel(model, options)
    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == :Optimal

    if sign(starting_point) < 0.0
        @test abs(JuMP.value.(model[:x]) - 1.0) < 1e-4
    else
        @test abs(JuMP.value.(model[:x]) + 1.0) < 1e-4
    end
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
function BatchLSModel(arg0::Vector{IntegratedPropagatorBuilder}, arg1::List, arg2::ParameterDriversList, arg3::ModelObserver)
    return BatchLSModel((Vector{IntegratedPropagatorBuilder}, List, ParameterDriversList, ModelObserver), arg0, arg1, arg2, arg3)
end

function create_propagators(obj::BatchLSModel, arg0::RealVector)
    return jcall(obj, "createPropagators", Vector{NumericalPropagator}, (RealVector,), arg0)
end

function fetch_evaluated_measurement(obj::BatchLSModel, arg0::jint, arg1::EstimatedMeasurement)
    return jcall(obj, "fetchEvaluatedMeasurement", void, (jint, EstimatedMeasurement), arg0, arg1)
end

function get_evaluations_count(obj::BatchLSModel)
    return jcall(obj, "getEvaluationsCount", jint, ())
end

function get_iterations_count(obj::BatchLSModel)
    return jcall(obj, "getIterationsCount", jint, ())
end

function get_selected_propagation_drivers_for_builder(obj::BatchLSModel, arg0::jint)
    return jcall(obj, "getSelectedPropagationDriversForBuilder", ParameterDriversList, (jint,), arg0)
end

function is_forward_propagation(obj::BatchLSModel)
    return jcall(obj, "isForwardPropagation", jboolean, ())
end

function set_evaluations_counter(obj::BatchLSModel, arg0::Incrementor)
    return jcall(obj, "setEvaluationsCounter", void, (Incrementor,), arg0)
end

function set_iterations_counter(obj::BatchLSModel, arg0::Incrementor)
    return jcall(obj, "setIterationsCounter", void, (Incrementor,), arg0)
end

function value(obj::BatchLSModel, arg0::RealVector)
    return jcall(obj, "value", Pair, (RealVector,), arg0)
end

ext: jl


text:
# This file illustrates how the delta method may be used
# to calculate the covariance of a nonlinear function of
# estimated parameters
using Calculus, LinearAlgebra, Statistics
# this function that gives the elasticities of x*b wrt x:
function ElasticitiesLinearModel(theta, x)
	elasticities = (theta .* x) / (x'theta)
end

function main()
# define some random data
n = 100
k = 5
x = [ones(n,1) rand(n,k-1)]
beta = (1:k) # true betas are [1,2,3,4,5]
y = x*beta + 3.0*randn(n)
b, varb, junk = ols(y,x)

# we'll evaluate elasticities at the sample means of the data
x = vec(mean(x,dims=1))

Elasticities = theta -> vec(ElasticitiesLinearModel(theta, x))
# want derivatives by param in the rows: kXg from GMM theory
J = Calculus.jacobian(Elasticities, b, :central)'
elasticities = Elasticities(b)
var_elasticities = J*varb*J'
# Taa daa - the results!
println()
println("Delta method elasticities and standard errors, at the means of the regressors")
prettyprint([elasticities diag(var_elasticities)])
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# ---
# title: 726. Number of Atoms
# id: problem726
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Stack, Recursion
# link: <https://leetcode.com/problems/number-of-atoms/description/>
# hidden: true
# ---
#
# Given a chemical `formula` (given as a string), return the count of each atom.
#
# The atomic element always starts with an uppercase character, then zero or
# more lowercase letters, representing the name.
#
# One or more digits representing that element's count may follow if the count
# is greater than 1. If the count is 1, no digits will follow. For example, H2O
# and H2O2 are possible, but H1O2 is impossible.
#
# Two formulas concatenated together to produce another formula. For example,
# H2O2He3Mg4 is also a formula.
#
# A formula placed in parentheses, and a count (optionally added) is also a
# formula. For example, (H2O2) and (H2O2)3 are formulas.
#
# Given a `formula`, return _the count of all elements as a string in the
# following form_ : the first name (in sorted order), followed by its count (if
# that count is more than 1), followed by the second name (in sorted order),
# followed by its count (if that count is more than 1), and so on.
#
#
#
#
#
# **Example 1:**
#
#
#
#     Input: formula = "H2O"
#     Output: "H2O"
#     Explanation: The count of elements are {'H': 2, 'O': 1}.
#
#
# **Example 2:**
#
#
#
#     Input: formula = "Mg(OH)2"
#     Output: "H2MgO2"
#     Explanation: The count of elements are {'H': 2, 'Mg': 1, 'O': 2}.
#
#
# **Example 3:**
#
#
#
#     Input: formula = "K4(ON(SO3)2)2"
#     Output: "K4N2O14S4"
#     Explanation: The count of elements are {'K': 4, 'N': 2, 'O': 14, 'S': 4}.
#
#
# **Example 4:**
#
#
#
#     Input: formula = "Be32"
#     Output: "Be32"
#
#
#
#
# **Constraints:**
#
#   * `1 <= formula.length <= 1000`
#   * `formula` consists of English letters, digits, `'('`, and `')'`.
#   * `formula` is always valid.
#
#
## @lc code=start
using LeetCode

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


text:
function resolve_colnames(l, r, on)
    isnothing(on) && ArgumentError("Missing join argument 'on'")
    if on isa Symbol
        l_cmp_symbols = [on]
        r_cmp_symbols = [on]
    elseif on isa Pair{Symbol,Symbol}
        l_cmp_symbols = getindex(on, 1)
        r_cmp_symbols = getindex(on, 2)
    elseif on isa Vector{Symbol}
        l_cmp_symbols = on
        r_cmp_symbols = on
    elseif on isa Vector{Pair{Symbol,Symbol}}
        l_cmp_symbols = getindex.(on, 1)
        r_cmp_symbols = getindex.(on, 2)
    end

    l_cols = Tables.columns(l)
    l_schema = Tables.schema(l_cols)
    cmp_l_indices = Tables.columnindex.(Ref(l_schema), l_cmp_symbols)
    other_l_indices = setdiff(1:length(l_cols), cmp_l_indices)
    cmp_l = NTuple{length(cmp_l_indices),Int}(cmp_l_indices)
    other_l = NTuple{length(other_l_indices),Int}(other_l_indices)

    r_cols = Tables.columns(r)
    r_schema = Tables.schema(r_cols)
    cmp_r_indices = Tables.columnindex.(Ref(r_schema), r_cmp_symbols)
    other_r_indices = setdiff(1:length(r_cols), cmp_r_indices)
    cmp_r = NTuple{length(cmp_r_indices),Int}(cmp_r_indices)
    other_r = NTuple{length(other_r_indices),Int}(other_r_indices)

    rnames = (name in l_schema.names ? Symbol(string(name) * "_2") : name for name in getindex.(Ref(r_schema.names), other_r_indices))
    final_nameset = (l_schema.names..., rnames...)

    final_nameset, other_l, other_r, cmp_l, cmp_r
end

"""
    match_inner_indices(l, r, l_ind::NTuple{N,Int}, r_ind::NTuple{N,Int})

Returns two vectors containing indices of matched rows.
Standard non-optimized use case.
"""
function match_inner_indices(l, r, l_ind::NTuple{N,Int}, r_ind::NTuple{N,Int}) where {N}
    l_length = length(Tables.rows(l))
    vl = Vector{UInt}()
    vr = Vector{UInt}()
    sizehint!(vl, l_length)
    sizehint!(vr, l_length)
    for (oind, oel) in enumerate(Tables.rows(l))
        for (iind, iel) in enumerate(Tables.rows(r))
            if compare_rows_eq(oel, iel, l_ind, r_ind)
                push!(vl, oind)
                push!(vr, iind)
            end
        end
    end
    vl, vr
end

"""
    match_inner_indices_lookup(l, lookup, l_ind::NTuple{N,Int})

Returns two vectors containing indices of matched rows.
Uses `lookup` to find the matching indices.

`lookup` needs to be a dict-like structure that contains keys in
form of a `Tuple` of all matching columns and values in form
of type `Vector{UInt}` containing the related row indices.
"""
function match_inner_indices_lookup(l, lookup, l_ind::NTuple{N,Int}) where {N}
    l_length = length(Tables.rows(l))
    vl = Vector{UInt}()
    vr = Vector{UInt}()
    sizehint!(vl, l_length)
    sizehint!(vr, l_length)

    row_tuple = (row, cols) -> ([Tables.getcolumn(row, x) for x in cols]...,)
    _evec = Vector{UInt}()
    for (oind, oel) in enumerate(Tables.rows(l))
        indices = get(lookup, row_tuple(oel, l_ind), _evec)
        for iind in indices
            push!(vl, oind)
            push!(vr, iind)
        end
    end
    vl, vr
end

"""
    match_inner_indices_lsorted_rsorted(l, r, cmp_l::NTuple{N,Int}, cmp_r::NTuple{N,Int}, runique::Bool)

Returns two vectors containing indices of matched rows.
Optimized pass for the left table sorted, right table sorted and optionally right table only containing unique keys.
"""
function match_inner_indices_lsorted_rsorted(l, r, cmp_l::NTuple{N,Int}, cmp_r::NTuple{N,Int}, runique::Bool) where {N}
    l_length = length(Tables.rows(l))
    vl = Vector{UInt}()
    vr = Vector{UInt}()
    sizehint!(vl, l_length)
    sizehint!(vr, l_length)

    ri = enumerate(Tables.rows(r))
    li = enumerate(Tables.rows(l))

    riter = iterate(ri)
    liter = iterate(li)
    while true
        if riter === nothing || liter === nothing
            break
        end

        (iind, iel), rstate = riter
        (oind, oel), lstate = liter

        if compare_rows_lt(iel, oel, cmp_r, cmp_l)
            riter = iterate(ri, rstate)
            continue
        end

        if compare_rows_lt(oel, iel, cmp_l, cmp_r)
            liter = iterate(li, lstate)
            continue
        end

        if compare_rows_eq(oel, iel, cmp_l, cmp_r)
            push!(vl, oind)
            push!(vr, iind)
            if runique
                liter = iterate(li, lstate)
                continue
            else
                inner_riter = iterate(ri, rstate)
                while inner_riter !== nothing
                    (i_iind, i_iel), i_rstate = inner_riter
                    if compare_rows_eq(oel, i_iel, cmp_l, cmp_r)
                        push!(vl, oind)
                        push!(vr, i_iind)
                        inner_riter = iterate(ri, i_rstate)
                    else
                        break
                    end
                end
                liter = iterate(li, lstate)
            end
        end
    end
    vl, vr
end

"""
    match_inner_indices_runique(l, r, cmp_l::NTuple{N,Int}, cmp_r::NTuple{N,Int})

Returns two vectors containing indices of matched rows.
Optimized pass for joins with the right table containing unique keys only.
"""
function match_inner_indices_runique(l, r, cmp_l::NTuple{N,Int}, cmp_r::NTuple{N,Int}) where {N}
    l_length = length(Tables.rows(l))
    vl = Vector{UInt}()
    vr = Vector{UInt}()
    sizehint!(vl, l_length)
    sizehint!(vr, l_length)

    for (oind, oel) in enumerate(Tables.rows(l))
        for (iind, iel) in enumerate(Tables.rows(r))
            if compare_rows_eq(oel, iel, cmp_l, cmp_r)
                push!(vl, oind)
                push!(vr, iind)
                break
            end
        end
    end
    vl, vr
end

"""
    match_inner_indices_rsorted(l, r, cmp_l::NTuple{N,Int}, cmp_r::NTuple{N,Int})

Returns two vectors containing indices of matched rows.
Optimized pass for joins with a sorted right table.
"""
function match_inner_indices_rsorted(l, r, cmp_l::NTuple{N,Int}, cmp_r::NTuple{N,Int}) where {N}
    l_length = length(Tables.rows(l))
    vl = Vector{UInt}()
    vr = Vector{UInt}()
    sizehint!(vl, l_length)
    sizehint!(vr, l_length)

    for (oind, oel) in enumerate(Tables.rows(l))
        prev_cmp = false
        for (iind, iel) in enumerate(Tables.rows(r))
            cmp = compare_rows_eq(oel, iel, cmp_l, cmp_r)
            if cmp
                push!(vl, oind)
                push!(vr, iind)
            end
            prev_cmp && !cmp && break
            prev_cmp = cmp
        end
    end
    vl, vr
end

"""
    find_outer_indices(d, inner_indices)

Finds the unmatched indices from the table.
"""
function find_outer_indices(d, inner_indices)
    s = Set(one(UInt):length(Tables.rows(d)))
    setdiff!(s, inner_indices)
end

"""
    build_joined_table(jointype, names, l, r, inner_l, inner_r, outer_l, other_r)

Takes the indices of matching rows (`inner*`) and the ones that weren't matched (`outer_l`) from the `l` table
and builds the result based on that.

Uses all the columns from the left column and the `other_r` columns from the right table.
"""
function build_joined_table(
        jointype::Symbol,
        names::Tuple,
        l,
        r,
        inner_l::Vector{UInt},
        inner_r::Vector{UInt},
        outer_l::Set{UInt},
        other_r::NTuple,
    )

    fulllength = length(inner_l) + length(outer_l)

    cols = Vector{AbstractVector}(undef, length(names))

    colcounter = one(Int)

    for c in Tables.columns(l)
        newc = Vector{eltype(c)}(undef, fulllength)
        copyto!(newc, view(c, inner_l))
        copyto!(newc, length(inner_l) + 1, view(c, collect(outer_l)))
        cols[colcounter] = newc
        colcounter += 1
    end

    rcols = Tables.columns(r)
    for i in other_r
        t = Tables.schema(rcols).types[i]
        vectype = jointype == :innerjoin ? t : Union{t,Missing}
        newc = Vector{vectype}(undef, fulllength)
        if length(inner_r) > 0 # skip fetching and copying if there's no records matched
            c = Tables.getcolumn(rcols, i)
            copyto!(newc, view(c, inner_r))
        end
        cols[colcounter] = newc
        colcounter += 1
    end

    sink = Tables.materializer(l)
    sink((; zip(names, cols)...))
end

@inline function compare_rows_eq(o, i, cmp_l::NTuple{N,Int}, cmp_r::NTuple{N,Int}) where {N}
    test = true
    @inbounds for x = 1:N
        test &= Tables.getcolumn(o, cmp_l[x]) == Tables.getcolumn(i, cmp_r[x])
        test || break
    end
    test
end

@inline function compare_rows_lt(o, i, cmp_l::NTuple{N,Int}, cmp_r::NTuple{N,Int}) where {N}
    @inbounds for x = 1:N
        l = Tables.getcolumn(o, cmp_l[x])
        r = Tables.getcolumn(i, cmp_r[x])
        if l > r
            return false
        elseif l < r
            return true
        end
    end
    false
end


text:
using ConstraintSolver
CS = ConstraintSolver

function sudokus_from_file(filename, sep='\n')
    s = open(filename) do file
        read(file, String)
    end
    str_sudokus = split(strip(s), sep)
    grids = AbstractArray[]
    for str_sudoku in str_sudokus
        str_sudoku = strip(replace(str_sudoku, "."=>"0"))
        if length(str_sudoku) != 81
            continue
        end
        one_line_grid = parse.(Int, split(str_sudoku,""))
        grid = reshape(one_line_grid, 9, 9)
        push!(grids, grid)
    end
    return grids
end

function create_sudoku_grid!(com, grid)
    com_grid = Array{CS.Variable, 2}(undef, 9, 9)
    for (ind,val) in enumerate(grid)
        if val == 0
            com_grid[ind] = add_var!(com, 1, 9)
        else
            com_grid[ind] = add_var!(com, 1, 9; fix=val)
        end
    end
    return com_grid
end

function add_sudoku_constr!(com, grid)
    for rc=1:9
        #row
        variables = grid[CartesianIndices((rc:rc,1:9))]
        add_constraint!(com, CS.all_different([variables...]))
        #col
        variables = grid[CartesianIndices((1:9,rc:rc))]
        add_constraint!(com, CS.all_different([variables...]))
    end

    for br=0:2
        for bc=0:2
            variables = grid[CartesianIndices((br*3+1:(br+1)*3,bc*3+1:(bc+1)*3))]
            add_constraint!(com, CS.all_different([variables...]))
        end
    end
end

function fulfills_sudoku_constr(com_grid)
    correct = true
    for rc=1:9
        row = com_grid[CartesianIndices((rc:rc,1:9))]
        if any(v->!CS.isfixed(v), row)
            return false
        end
        vals = unique([value(v) for v in row])
        correct = length(vals) != 9 ? false : correct

        col = com_grid[CartesianIndices((1:9, rc:rc))]
        if any(v->!CS.isfixed(v), col)
            return false
        end
        vals = unique([value(v) for v in col])
        correct = length(vals) != 9 ? false : correct
    end

    for br=0:2
        for bc=0:2
            box = com_grid[CartesianIndices((br*3+1:(br+1)*3,bc*3+1:(bc+1)*3))]
            if any(v->!CS.isfixed(v), box)
                return false
            end
            vals = unique([value(v) for v in box])
            correct = length(vals) != 9 ? false : correct
        end
    end
    return correct
end

function print_search_space(com_grid; max_length=:default)
    if max_length == :default
        if all(v->CS.isfixed(v), com_grid)
            max_length = 2
        else
            max_length = 20
        end
    end

    for y=1:size(com_grid)[1]
        line = ""
        for x=1:size(com_grid)[2]
            if !CS.isfixed(com_grid[y,x])
                possible = sort(values(com_grid[y,x]))
                pstr = join(possible, ",")
                space_left  = floor(Int, (max_length-length(pstr))/2)
                space_right = ceil(Int, (max_length-length(pstr))/2)
                line *= repeat(" ", space_left)*pstr*repeat(" ", space_right)
            else
                pstr = string(value(com_grid[y,x]))
                space_left  = floor(Int, (max_length-length(pstr))/2)
                space_right = ceil(Int, (max_length-length(pstr))/2)
                line *= repeat(" ", space_left)*pstr*repeat(" ", space_right)
            end
        end
        println(line)
    end
end


text:
module BOMCP

using Plots

using Random
using Statistics
using StatsBase
using LinearAlgebra
using NLopt
using CPUTime
using Printf
using ElasticArrays
using PDMats
using NearestNeighbors
using Distributions
using Kalman
import PDMats: *, \, diag

using GaussianProcesses
import GaussianProcesses: MeanConst, predict_f
using BayesianOptimization
import BayesianOptimization: AbstractAcquisition

using POMDPs
using POMDPPolicies
using POMDPModelTools
using POMDPSimulators
using BeliefUpdaters

export
    BOMCPTree,
    insert_state_node!,
    insert_action_node!,
    BOMCPStateNode,
    children,
    n_children,
    isroots
include("./tree.jl")

export
    estimate_value,
    RolloutEstimator,
    SolvedRolloutEstimator,
    convert_estimator,
    convert_to_policy,
    rollout,
    init_Q,
    init_N
include("./node_initialization.jl")

export
    GaussianObsUpdater,
    EKFUpdater,
    x2s, s2x, gen_A, gen_C
include("./beliefs.jl")

export
    ReportWhenUsed,
    default_action
include("./default_action.jl")

export
    GPLA,
    reset_gp!,
    predict_fix,
    predict_test
include("./gaussian_process.jl")

export
    rand_action,
    build_gp,
    ActionSelector,
    BOActionSelector
include("./action_selector.jl")

export
    update_gp!,
    next_action,
    next_ow_action,
    vectorize,
    build_action_set
include("./bayes_opt.jl")

export
    BOMCPSolver,
    BOMCPPlanner
include("./solver.jl")

export
    clear_tree!,
    simulate,
    action,
    action_info
include("./search.jl")

include("./visualization.jl")

end


text:
type bracket{n}
    i::Vector{Int}
    w::Float64
end

function bracket(v::Float64, x::Vector)
    dx = (x[end] - x[1]) / (length(x) - 1)
    i,d = divrem(v - x[1], dx)
    if d == 0
        return [bracket{1}([Int(i + 1)], (1 - d / dx))]
    elseif d == 1
        return [bracket{1}([Int(i + 2)], d / dx)]
    else
        return [bracket{1}([Int(i + 1)], (1 - d / dx)),bracket{1}([Int(i + 2)], d / dx)]
    end
end

promote_type{n1, n2}(::Type{bracket{n1}}, ::Type{bracket{n2}}) = Type(bracket{n1 + n2})
*{n1, n2}(a::bracket{n1}, b::bracket{n2}) = bracket{n1 + n2}(vcat(a.i, b.i), a.w * b.w)
*{n2}(a::bracket{0}, b::bracket{n2}) = b
kron{n}(a::Array{bracket{n}, 1}) = a

convert(::Type{Tuple}, x::Vector{Int}) = ntuple(i -> x[i], length(x))
spzeros(n::Int) = SparseVector(n, Int[], Float64[])

"
   updateT(M)

Compute transition matrix for model M using tensor grid.
"
function updateT(M::Model)
    N = length(M.distribution.G[1])
    Endo = Endogenous(M.variables)
    Stoc = Stochastic(M.variables)
    Pf = Vector{Float64}[clamp!(M(s.name, hcat([g[:] for g in M.distribution.G]...)), s.bounds...) for s in Endogenous(M.variables)]
    Gdims = size(M.distribution.G[1])
    Sdim   = tuple([length(s) for s in Stochastic(M.variables)]...)

    ns = length(Stoc)
    nT = prod(Sdim) * length(Endo) * 2
    vals = zeros(N, nT)
    rloc = zeros(UInt32, N, nT)
    cloc = zeros(UInt32, N, nT)
    ccnt = zeros(UInt32, N)

    ne = length(Endo)

    @threads for i = 1:N
        brackets = [bracket(Pf[ie][i], Endo[ie].x) for ie = 1:length(Endo)]
        S1 = collect(ind2sub(Gdims, i))[ne + 1:end]
        abrackets = kron(brackets...)

        cnt = 0
        for j = 1:length(abrackets)
            for jj ∈ CartesianRange(Sdim)
                t = 1.0
                for S2 = 1:ns
                    if isag(Stoc[S2])
                        t *= Stoc[S2].T[jj.I[S2], S1[S2]]
                    else
                        t *= Stoc[S2].T[S1[S2], jj.I[S2]]
                    end
                end
                cnt += 1
                rloc[i,cnt] = sub2ind(Gdims, abrackets[j].i..., jj.I...)
                vals[i,cnt] = abrackets[j].w * t
                cloc[i,cnt] = i
                ccnt[i] = cnt
            end
        end
    end
    rloc, cloc, vals = vec(rloc), vec(cloc), vec(vals)
    id = (vals .== 0)
    M.distribution.T = sparse(rloc[!id], cloc[!id], vals[!id], N, N)
    return
end

function updated(M::Model)
    d0 = M.distribution.d[:] + .25
    for j = 1:5000
        d1 = M.distribution.T * d0
        mean(abs.((d1[:] - d0[:]))) < 1e-10 && break
        d0[:] = d1[:] / sum(d1)
        if j == 5000
            warn("Model distribution did not converge")
        end
    end

    M.distribution.d[:] = d0
    return
end

function updatedistribution(M::Model)
    updateT(M)
    updated(M)
end

function ∫(M::Model, v::Symbol, t = 0)
    isa(M.variables[v], Dependant) && M.variables[v].update(M)
    ∫(M, M[v, t])
end

function ∫(M::Model, v::Vector{Float64})
    ∫(M, reshape(M.G(v, hcat([vec(x) for x in M.distribution.G]...)), size(M.distribution.G[1])))
end

function ∫(M::Model, V::Array{Float64})
    d = V .* M.distribution.d
    if all(!Bool[isag(s) for s in State(M.variables)])
        return sum(d)
    else
        ns = length(State(M.variables))
        id = find([isag(s) for s in State(M.variables)])
        D = zeros(size(M.distribution.d, id...))
        for i = 1:prod(size(D))
            id1 = ind2sub(size(D), i)
            inds = Any[Colon() for i = 1:ns]
            cnt = 1
            for ii = 1:ns
                if [isag(s) for s in State(M.variables)][ii]
                    inds[ii] = id1[cnt]
                    cnt += 1
                end
            end

            D[id1...] = sum(d[inds...]) / max(sum(M.distribution.d[inds...]), 1e-16)
        end
        return D
    end
end

function updateAggregates(M::Model, ϕ = 1.0)
    if any(Bool[isag(s) for s in State(M.variables)])
        agg = Aggregate(M.variables)
        for s ∈ Aggregate(M.variables)
            ag  = ∫(M, M[s.target.name, timeof(s.target)])
            eag = filter(x -> isag(x), State(M.variables))
            t = [isa(e, Endogenous) ? -1 : 0 for e in eag]

            X   = ndgrid([e.x for e in eag]...)
            Xid = findfirst(M.variables, s)
            # tX = zeros(size(M.X,1))
            for ii = 1:length(X[1])
                id = BitArray(.*([M[eag[ie].name,t[ie]] .== X[ie][ii] for ie = 1:length(eag)]...))
                M.X[id,Xid] *= (1 - ϕ)
                M.X[id,Xid] += ϕ * ag[ii] * ones(sum(id))
            end
        end
    else
        for s ∈ Aggregate(M.variables)
            # M.X[:,findfirst(M.variables,s)] *= (1-ϕ)
            # M.X[:,findfirst(M.variables,s)] += ϕ*∫(M,s.target.name,timeof(s.target))
            M.X[:,findfirst(M.variables, s)] = (1 - ϕ) * M.X[:, findfirst(M.variables, s)] + ϕ * ∫(M, s.target.name, timeof(s.target))
        end
    end
end

function setaggregate!(M::Model, s::Symbol, x)
    if all(!Bool[isag(s) for s in State(M.variables)])
        M[s,0] = x
    else
        eag = filter(x -> isag(x), State(M.variables))
        t = [isa(e, Endogenous) ? -1 : 0 for e in eag]
        X   = ndgrid([e.x for e in eag]...)
        @assert size(X[1]) == size(x) "size(x) must equal size of aggregate state space"
        # Xid = findfirst(M.variables,M[s])
        Xid = findfirst(x -> x.name == s && !isa(x, State), M.variables)
        for ii = 1:length(X[1])
            id = BitArray(.*([M[eag[ie].name,t[ie]] .== X[ie][ii] for ie = 1:length(eag)]...))
            M.X[id,Xid] = x[ii]
        end
    end
end

function updateA(M, ϕ = 1.0)
    updateallD(M)
    updatedistribution(M)
    updateAggregates(M, ϕ)
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
module Mixins

export @use_mixins, @define_mixin

using MacroTools: prewalk, @capture
using MLStyle: @match

macro use_mixins(code)
  mixins = Symbol[]
  stripped_code = prewalk(code) do exp
    @match exp begin
      Expr(args...) => begin
        for arg in args
          @match arg begin
            Expr(:macrocall, macroargs...) => begin
              if macroargs[1] == Symbol("@mixin")
                push!(mixins, macroargs[3])
              end
            end
            _ => ()
          end
        end
        exp
      end
      _ => exp
    end
  end
  unique!(mixins)
  @match mixins begin
    [] => esc(stripped_code)
    [m,rest...] => esc(Expr(:macrocall, Symbol("@" * string(m)), nothing, [], copy(rest), stripped_code))
  end
end

macro expand_mixins(mixin_defs, code)
  mixin_defs = Dict(mixin_defs...)
  unescaped = prewalk(code) do exp
    @match exp begin
      Expr(head, args...) => begin
        newargs = map(args) do arg
          @match arg begin
            Expr(:macrocall, macroargs...) => begin
              if macroargs[1] == Symbol("@mixin")
                mixin_defs[macroargs[3]]
              end
            end
            _ => [arg]
          end
        end
        Expr(head, vcat(newargs...)...)
      end
      _ => exp
    end
  end
  esc(unescaped)
end

function shift_mixins(newdef,defs,remaining,code)
  newdefs = [defs...;newdef]
  if length(remaining) == 0
    Expr(:macrocall, GlobalRef(Mixins, Symbol("@expand_mixins")), nothing, newdefs, code)
  else
    m,rest... = remaining
    Expr(:macrocall, Symbol("@"*string(m)), nothing, newdefs, copy(rest), code)
  end |> esc
end

macro define_mixin(name, block)
  newdef = name => copy(block.args)
  :(macro $(esc(name))(defs, remaining, code)
      $(GlobalRef(Mixins,:shift_mixins))($newdef, defs, remaining, code)
    end)
end

end # module


text:
module SphericalFunctions

using FastTransforms, LinearAlgebra, ProgressMeter, Quaternionic
import Base.Threads: @threads, nthreads

const MachineFloat = Union{Float16, Float32, Float64}

include("utils.jl")

include("complex_powers.jl")
export complex_powers, complex_powers!

include("indexing.jl")
export Ysize, Yrange, Yindex, deduce_limits, theta_phi, phi_theta
export WignerHsize, WignerHindex, _WignerHindex, WignerHrange
export WignerDsize, WignerDindex, WignerDrange

include("iterators.jl")
export Diterator, diterator, Yiterator

include("associated_legendre.jl")
export ALFRecursionCoefficients, ALFrecurse!, ALFcompute!, ALFcompute

include("Hrecursion.jl")
export H!, H_recursion_coefficients

include("evaluate.jl")
export d!, d, D!, Y!
export dprep, dstorage, Dprep, Dstorage, Yprep, Ystorage

include("weights.jl")
export fejer1, fejer2, clenshaw_curtis

include("map2salm.jl")
export map2salm, map2salm!, plan_map2salm

#include("rotate.jl")
#export rotate!

end # module


text:
# ====================  Auxiliary functions ====================

"""
    build_adj_matrix(pred, indl)

Constructs adjacency matrix with vertices from `indl` and edges
from `pred(i,j)` where `i,j ∈ indl`
"""
function build_adj_matrix(pred::Predicate, indl)
  adj = BitArray(undef, length(indl),length(indl))
  for (i,ind) in enumerate(indl)
    adj[i,i] = true
    for (j,ind2) in enumerate(Iterators.drop(indl,i))
      adj[i,j+i] = pred(ind,ind2) || pred(ind2,ind)
      adj[j+i,i] = adj[i,j+i]
    end
  end
  return adj
end

"""
    find_classes(adj::BitArray{2})

returns an array of length `size(adj,1)` with each entry with index `i` being
a unique identifier for the equivalency class of the node `i`.
Uses adjacency matrix as input. I.e. graph needs to be constructed before.

# Examples
```
julia> find_classes(convert(BitArray, [0 1 0; 1 0 0; 0 0 0]))
[1, 1, 2]
```
"""
function find_classes(adj::BitArray{2})
  classes = collect(1:size(adj,1))              # every node has its own class at the start
  cc = 1                                        # current class counter
  openL = trues(size(adj,1))                    # no vertices have been visited yet
  for i in 1:size(adj,1)                        # check all vertices at least once
    !openL[i] && continue                       # if node already checked, continue
    classes[i] = cc                             # current node has class cc
    indL = [i]                                  # add current vertex to search list
    while !isempty(indL)                        # search through all reachable vertices
      j = pop!(indL)                            # first vertex in index List is next candidate
      !openL[j] && continue
      openL[j] = false                          # set current node to visited
      neighbors = findall(adj[i,:] .& openL)    # add all adjacent AND open vertices to class
      indL = union(neighbors,indL)              # without double entries
      classes[indL] .= cc                       # set all neighbors to current class
    end
    cc += 1
    sum(openL) == 0 && break                    # no more open vertices left, return
  end
  return classes
end

"""
    find_classes(m::Mapping, vl)

returns an array of length `size(adj,1)` with each entry with index `i` being
a unique identifier for the equivalency class of the node `i`.
Uses a mapping function `m(x) = y` and an vertex list `vl` to construct graph i
n place.
`m` has to be symmetric (i.e. `m(x) = y => m(y) = x`)!
In case the vertex list is closed under the mapping (i.e. there exist no `x` in
`vl` such that `m(x)` is not in `vl`), `closed` can be set to `false` in order
to improve performace.

# Examples
```
julia> find_classes(Mapping(x-> [-x]), -2:4)
Dict(0 => 3,4 => 5,2 => 1,-2 => 1,-1 => 2,3 => 4,1 => 2)
```
"""
function find_classes(m::Mapping, vl::AbstractArray; closed=false)
  classes = Dict(zip(vl,1:length(vl)))
  cc = 1                                        # current class counter
  openL = Dict(zip(vl,trues(length(vl))))       # mark all vertices as open (not visited)
  for vi in vl                                  # visit each entry at least once
    !openL[vi] && continue                      # if entry already checked, continue
    classes[vi] = cc                            # current node has class cc
    searchL = [vi]                              # add current vertex to search list
    while !isempty(searchL)                     # search through all reachable entries
      vj = pop!(searchL)                        # first vertex in index List is next candidate
      !openL[vj] && continue
      openL[vj] = false                         # set current node to visited
      neighbors = if closed
        filter(x-> openL[x], m(vj))    # add all adjacent AND open vertices to class
      else
        filter(x->(x in keys(openL)) && openL[x], m(vj))    # add all adjacent AND open vertices to class, skip entries outside vl
      end
      searchL = union(neighbors,searchL)        # without double entries
      for el in searchL
          classes[el] = cc
      end
    end
    cc += 1
    sum(values(openL)) == 0 && break                    # no more open vertices left, return
  end
  return classes
end

"""
    find_classes(pred::Predicate, indl; isSymmetric=false)

returns an array of length `length(indl)` with each entry with index `i` being
a unique identifier for the equivalency class of the node `i`.
"""
function find_classes(pred::Predicate, indl::IndexList; isSymmetric=false)
  adj = isSymmetric ? build_adj_matrix(pred, indl) :
    build_adj_matrix(Predicate((x,y) ->(pred(x,y) || pred(y,x))), indl)
  return find_classes(adj)
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie
GLMakie.activate!() # HIDE
let
    u = LinRange(-1.5, 2, 50)
    v = LinRange(0, 2*pi, 50)
    X1 = [u for u in u, v in v]
    Y1 = [(u^2 + 1)*cos(v) for u in u, v in v]
    Z1 = [(u^2 + 1)*sin(v)  for u in u, v in v]
    fig, _ = surface(X1, Y1, Z1, shading = true, ambient = Vec3f0(0.95, 0.95, 0.95),
        backlight = 0.95f0, color = sqrt.(X1.^2 .+ Y1.^2 .+ Z1.^2), colormap = :viridis)
    wireframe!(X1, Y1, Z1)
    fig
    save(joinpath(@__DIR__, "output", "revolutionSurf2.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status("GLMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
"""
Represents any kind of wrapper structure that was generated from a Vulkan structure.
`D` is a `Bool` parameter indicating whether the structure has specific dependencies or not.
"""
abstract type VulkanStruct{D} end

Base.broadcastable(x::VulkanStruct) = Ref(x)

const FunctionPtr = Union{Ptr{Cvoid}, Base.CFunction}
const OptionalPtr{T} = Union{T, Ptr{Cvoid}}

"""
High-level structure with idiomatic Julia types.
"""
abstract type HighLevelStruct end

Base.broadcastable(x::HighLevelStruct) = Ref(x)

hl_type(T::Type{<:HighLevelStruct}) = T
intermediate_type(T::Type{<:VulkanStruct}) = T

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
module FieldApplyBenchs

using InplaceArrays.Arrays
using InplaceArrays.Fields
using InplaceArrays.Fields: MockField, MockBasis

@inline function repeat(n,f,args...)
  for i in 1:n
    f(args...)
  end
  nothing
end

function bench1(n)
  np = 4
  p = Point(1,2)
  x = fill(p,np)
  v = 3.0
  d = 2
  f = MockField{d}(v)
  g = apply_kernel_to_field(bcast(+),f,v)
  cg = field_cache(g,x)
  @time repeat(n,evaluate!,cg,g,x)
  ∇g = gradient(g)
  ∇cg = field_cache(∇g,x)
  @time repeat(n,evaluate!,∇cg,∇g,x)
  h = apply_kernel_to_field(bcast(+),f,f)
  ch = field_cache(h,x)
  @time repeat(n,evaluate!,ch,h,x)
end

#using InplaceArrays.Inference
#
#function bench3(n)
#  np = 4
#  p = Point(1,2)
#  x = fill(p,np)
#  v = 3.0
#  d = 2
#  f = MockField{d}(v)
#  k = bcast(+)
#  @time repeat(n,apply_kernel_to_field,k,f,v)
#end
#
#using InplaceArrays.Fields: Valued
#
#function bench4(n)
#  np = 4
#  p = Point(1,2)
#  x = fill(p,np)
#  v = 3.0
#  ndof = 8
#  w = fill(v,ndof)
#  d = 2
#  f = MockBasis{d}(v,ndof)
#  k = bcast(+)
#  fx = evaluate(f,x)
#  T = kernel_return_type(k,fx,fx)
#  @time repeat(n,apply_kernel_to_field,T,k,f,f)
#  @time repeat(n,apply_kernel_to_field,T,k,f,w)
#  @time repeat(n,apply_kernel_to_field,T,k,w,f)
#  vk = Valued(k,f,w)
#  cvk = kernel_cache(vk,f,w)
#  @time repeat(n,apply_kernel!,cvk,vk,f,w)
#end

for n in (1,1,10,1000,100000)
  @eval begin
    println("+++ runing suite for n = $($n) +++")
    bench1($n)
    #bench3($n)
    #bench4($n)
  end
end

end # module


text:
#Default options
const defaultParams =
    @compat Dict{Symbol,Any}(:plotlib => nothing,
                            :plotlib_set => false,
                            :storeresults => false,
                            :doc_number => 0,
                            :chunk_defaults =>
                                Dict{Symbol,Any}(
                                :echo=> true,
                                :results=> "markup",
                                :hold => false,
                                :fig=> true,
                                :include=> true,
                                :eval => true,
                                :tangle => true,
                                :cache => false,
                                :fig_cap=> nothing,
                                #Size in inches
                                :fig_width => 6,
                                :fig_height => 4,
                                :fig_path=> "figures",
                                :dpi => 96,
                                :term=> false,
                                :display => false,
                                :prompt => "\njulia> ",
                                :name=> nothing,
                                :wrap=> true,
                                :line_width => 75,
                                :engine=> "julia",
                                #:option_AbstractString=> "",
                                #Defined in formats
                                :fig_ext => nothing,
                                :fig_pos=> nothing,
                                :fig_env=> nothing,
                                :out_width=> nothing,
                                :out_height=> nothing,
                                )
                            )
#This one can be changed at runtime, initially a copy of defaults
const rcParams = deepcopy(defaultParams)

#Parameters set per document
const docParams =Dict{Symbol,Any}(
                                :fig_path=> nothing,
                                :fig_ext => nothing,
                            )

"""
`set_chunk_defaults(opts::Dict{Symbol, Any})`

Set default options for code chunks, use get_chunk_defaults
to see the current values.

e.g. set default dpi to 200 and fig_width to 8

```
julia> set_chunk_defaults(Dict{Symbol, Any}(:dpi => 200, fig_width => 8))
```
"""
function set_chunk_defaults(opts::Dict{Symbol, Any})
  merge!(rcParams[:chunk_defaults], opts)
  return nothing
end

"""
`get_chunk_defaults()`

Get default options used for code chunks.
"""
function get_chunk_defaults()
  return(rcParams[:chunk_defaults])
end

"""
`restore_chunk_defaults()`

Restore Weave.jl default chunk options
"""
function restore_chunk_defaults()
  rcParams[:chunk_defaults] = defaultParams[:chunk_defaults]
  merge!(rcParams[:chunk_defaults], docParams)
  return nothing
end


text:
# ---
# cover: example_2_2.png
# id: dip_example_2_2
# ---

# 本例展示了图片在不同空间分辨率下的显示效果

using Images, DIPTutorials

# 读取示例图片
ori_img = load(DIP3e.get_filepath("Fig0220")) # 3692×2812 Gray image

# 由于dpi是一个跟打印有关的概念，这里以原始图片为标准大小进行模拟：
#
# 1. 对原始图片进行下采样
# 2. 将下采样的图片再上采样回原尺寸
#
# 简单起见，下采样与上采样使用了`imresize`函数，关于`imresize`的实现在[例2-4](@ref dip_example_2_4)中会介绍

## 设置DPI
img_dpis = [1250, 300, 150, 72]
img_sizes = [ceil.(Int, dpi/1250 .* size(ori_img)) for dpi in img_dpis]
## 下采样
downsampled_imgs = [imresize(ori_img, sz) for sz in img_sizes]
## 上采样
upsampled_imgs = [imresize(img, size(ori_img)) for img in downsampled_imgs]
nothing #hide #!jl

# 展示结果：第一行为下采样之后的结果，第二行为对应放大之后显示的结果。
# 从结果可以看到，在相同尺寸下显示/打印出来时，图片分辨率越低意味着dpi越低，锯齿状也越明显
vcat(mosaicview(downsampled_imgs, nrow=1),
     mosaicview(upsampled_imgs, nrow=1))

# save cover #src
save("example_2_2.png", imresize(ori_img, (220, 200))) #src


text:
# module for testing each localization modules

module TestLocalization
  using Test

  # target modules
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/particle/particle.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/random_sampling/anime_mcl_rand_samp.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/systematic_sampling/anime_mcl_sys_samp.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/kld_sampling/particle_num_kld.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/kld_sampling/particle_num_wh.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/kld_sampling/anime_kld_mcl.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/parameter_adjustment/motion_test_forward/motion_test_forward.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/parameter_adjustment/motion_test_forward_bias/motion_test_forward_bias.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/parameter_adjustment/motion_test_rot_bias/motion_test_rot_bias.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/parameter_adjustment/sensor_test_static/sensor_test_static.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/occlusion_free_mcl/occlusion_free_particle.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/particle_filter/occlusion_free_mcl/anime_occlusion_free_mcl.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/extended_kalman_filter/anime_ekf.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/global_localization/test_global_kf.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/global_localization/test_global_mcl.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/global_localization/test_global_kld_mcl.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/reset_process/random_reset_mcl/anime_random_reset_mcl.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/reset_process/sensor_reset_mcl/anime_sensor_reset_mcl.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/reset_process/adaptive_reset_mcl/anime_adaptive_reset_mcl.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/reset_process/expansion_reset_mcl/anime_expansion_reset_mcl.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/localization/reset_process/expansion_sensor_reset_mcl/anime_expansion_sensor_reset_mcl.jl"))

  function main()
    @testset "ParticleFilter" begin
      @testset "Particle" begin
        p = Particle([0.0, 0.0, 0.0], 0.0)
        @test p.pose[1] == 0.0
        @test p.pose[2] == 0.0
        @test p.pose[3] == 0.0
        @test p.weight == 0.0
      end
      @testset "RandomSampling" begin
        @test_nowarn AnimeMclRandSamp.main(is_test=true)
      end
      @testset "SystematicSampling" begin
        @test_nowarn AnimeMclSysSamp.main(is_test=true)
      end
      @testset "KldSampling" begin
        @testset "ParticleNumKld" begin
          @test_nowarn ParticleNumKld.main(true)
        end
        @testset "ParticleNumWh" begin
          @test ParticleNumWh.num(0.1, 0.01, 2) == 34.0
          @test ParticleNumWh.num_wh(0.1, 0.01, 2) == 33.0
          @test_nowarn ParticleNumWh.main()
        end
        @testset "KldMcl" begin
          @test_nowarn AnimeKldMcl.main(is_test=true)
        end
        @testset "OcclusionFree" begin
          @testset "OcclusionFreeParticle" begin
            op = OcclusionFreeParticle([0.0, 0.0, 0.0], 0.0)
            @test op.pose[1] == 0.0
            @test op.pose[2] == 0.0
            @test op.pose[3] == 0.0
            @test op.weight == 0.0
          end
          @testset "OcclusionFreeMcl" begin
            @test_nowarn AnimeOcclusionFreeMcl.main(is_test=true)
          end
        end
      end
      @testset "ParameterAdjustment" begin
        @testset "MotionTestForward" begin
          @test_nowarn MotionTestForward.main()
        end
        @testset "MotionTestForwardBias" begin
          @test_nowarn MotionTestForwardBias.main()
        end
        @testset "MotionTestRotBias" begin
          @test_nowarn MotionTestRotBias.main()
        end
        @testset "SensorTestStatic" begin
          @test_nowarn SensorTestStatic.main()
        end
      end
    end
    @testset "KalmanFilter" begin
      @testset "ExtendedKalmanFilter" begin
        @test_nowarn AnimeEkf.main(is_test=true)
      end
    end
    @testset "GlobalLocalization" begin
      @testset "GlobalKf" begin
        @test_nowarn TestGlobalKf.main(5)
      end
      @testset "GlobalMcl" begin
        @test_nowarn TestGlobalMcl.main(5)
      end
      @testset "GlobalKldMcl" begin
        @test_nowarn TestGlobalKldMcl.main(5)
      end
    end
    @testset "ResetProcess" begin
      @testset "RandomResetMcl" begin
        @test_nowarn AnimeRandomResetMcl.main(is_test=true)
      end
      @testset "SensorResetMcl" begin
        @test_nowarn AnimeSensorResetMcl.main(is_test=true)
      end
      @testset "AdaptiveResetMcl" begin
        @test_nowarn AnimeAdaptiveResetMcl.main(is_test=true)
      end
      @testset "ExpansionResetMcl" begin
        @test_nowarn AnimeExpansionResetMcl.main(is_test=true)
      end
      @testset "ExpansionSensorResetMcl" begin
        @test_nowarn AnimeExpansionSensorResetMcl.main(is_test=true)
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# Tests for distances
function t_lda()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)

Test.@test try
	A05 |> (A05 |> j4pr.lda(true) )
	+A05 |> (A05 |> j4pr.lda(true) )

	A05 |> (A05 |> j4pr.lda(false) )
	+A05 |> (A05 |> j4pr.lda(false) )

	A05 |> (A05 |> j4pr.ldasub(true) )
	+A05 |> (A05 |> j4pr.ldasub(true) )

	A05 |> (A05 |> j4pr.ldasub(false) )
	+A05 |> (A05 |> j4pr.ldasub(false) )
	true
catch
	false
end

end


text:
module DrsMathProgSolverInterface

using Logging
@Logging.configure(level=DEBUG)

importall MathProgBase.SolverInterface

export DrsMathProgModel,
    DrsMathProgSolver,
    loadproblem!,
    optimize!,
    status,
    getreducedcosts,
    getconstrduals,
    getobjval,
    getsolution

immutable DrsMathProgSolver <: AbstractMathProgSolver
    options
end
DrsMathProgSolver(; kwargs...) = DrsMathProgSolver(kwargs)

type DrsMathProgModel <: AbstractLinearQuadraticModel
    A               # constraint coefficients
    b               # RHS
    c               # objective coefficients
    basis           # basis variables
    nonbasis        # nonbasis variables
    status          # solution status
    objval          # objective value
    solution        # solution
    maxiter         # max iteration number
end
LinearQuadraticModel(s::DrsMathProgSolver) = DrsMathProgModel(; s.options...)

function setparameters!(m::Union{DrsMathProgSolver, DrsMathProgModel}; kwargs...)
    for (option, value) in kwargs
        if option == :MaxIter
            m.maxiter = value
        elseif option == :Silent && value == true
            Logging.configure(level=OFF)
        elseif option == :LogLevel
            Logging.configure(level=value)
        else
            warn("option $option is unsupported")
        end
    end
end

function DrsMathProgModel(A=0, b=0, c=0; kwargs...)
    m = DrsMathProgModel(A, b, c, 0, 0, 0, 0, 0, 20)
    setparameters!(m; kwargs...)
    m
end

function loadproblem!(m::DrsMathProgModel, A, l, u, c, lb, ub, sense)
    @debug("loadproblem!: A $A, l $l, u $u, c $c, lb $lb, ub $ub, sense $sense")
    m.A = A
    m.b = zeros(size(ub))
    m.c = c

    DrsTransformToStandardForm!(m, lb, ub, sense)

    m.A = [1 m.c'; [zeros(length(m.b)) m.A]]
    m.b = [0; m.b]

    DrsFindPotentialBasis!(m)
end

function DrsFindPotentialBasis!(m::DrsMathProgModel)
    r, c = size(m.A)
    m.basis = zeros(Int, r)
    for ic in 1:c
        column = m.A[:,ic]
        if countnz(column) == 1
            ir = findfirst(x -> x == 1, column)
            if ir != 0 && m.basis[ir] == 0
                # add the column if current row has not been selected
                m.basis[ir] = ic
            end
        end
    end
    m.nonbasis = setdiff(1:c, m.basis)
end

function DrsTransformToStandardForm!(m::DrsMathProgModel, lb, ub, sense)
    @assert length(lb) == length(ub) """the lengths of
    lower bounds ($(length(lb))) and
    upper bounds ($(length(ub))) are different"""

    r, c = size(m.A)

    # check if b is negative

    for i in 1:length(lb)
        if lb[i] == typemin(typeof(lb[i]))
            # <
            if ub[i] < 0
                m.A[i,:] = -m.A[i,:]
                lb[i], ub[i] = -ub[i], -lb[i]
            end
        elseif ub[i] == typemax(typeof(ub[i]))
            # >
            if lb[i] < 0
                m.A[i,:] = -m.A[i,:]
                lb[i], ub[i] = -ub[i], -lb[i]
            end
        else
            # =
            if lb[i] < 0
                m.A[i,:] = -m.A[i,:]
                lb[i], ub[i] = -ub[i], -lb[i]
            end
        end
    end

    # add variables
    surplus = []

    for i in 1:length(lb)
        if lb[i] == typemin(typeof(lb[i]))
            # <, add slack
            m.A = [m.A zeros(r, 1)]
            m.A[i, end] = 1
            m.b[i] = ub[i]
            m.c = [m.c; 0]
        elseif ub[i] == typemax(typeof(ub[i]))
            # >, add surplus
            m.A = [m.A zeros(r, 1)]
            m.A[i, end] = -1
            m.b[i] = lb[i]
            m.c = [m.c; 0]
            push!(surplus, i)
        end
    end

    # =, add artificial
    for i in 1:length(lb)
        if lb[i] != typemin(typeof(lb[i])) && ub[i] != typemax(typeof(ub[i]))
            m.A = [m.A zeros(r, 1)]
            m.A[i, end] = 1
            m.b[i] = lb[i]
            m.c = [m.c; 0]
        end
    end

    # add artificial for surplus
    for i in surplus
        m.A = [m.A zeros(r, 1)]
        m.A[i, end] = 1
        m.c = [m.c; 0]
    end

    if sense == :Max
        m.c = -m.c
    end
end

function optimize!(m::DrsMathProgModel)
    @debug("A $(m.A)")
    @debug("b $(m.b)")
    @debug("c $(m.c)")
    @debug("basis $(m.basis)")
    @debug("nonbasis $(m.nonbasis)")

    nonbasis_orig = m.nonbasis[:]

    B = m.A[:,m.basis]
    @debug("B $B")

    invB = inv(B)
    @debug("invB $invB")

    iter = 0
    while true
        @debug("ITERATION $iter")

        @debug("basis $(m.basis)")
        @debug("nonbasis $(m.nonbasis)")

        N = m.A[:,m.nonbasis]
        @debug("N $N")

        # B = m.A[:,m.basis]
        # @debug("B $B")

        #invB = inv(B)
        @debug("invB $invB")

        @debug("b $(m.b)")

        delta = invB[1,:]' * N
        @debug("delta $delta")

        incoming = indmin(delta)
        @debug("incoming $incoming")

        if delta[incoming] >= 0
            @debug("optimal")
            m.status = :Optimal
            break
        end

        Xk = invB * N[:,incoming]
        @debug("Xk $Xk")

        theta = m.b[2:end] ./ Xk[2:end]
        @debug("theta $theta")

        # only positive
        positive = [i for i in theta if i > 0]
        idx = indmin(positive)
        outgoing = findfirst(x -> x == positive[idx], theta) + 1
        @debug("outgoing $outgoing")

        rect = [invB[:,2:end] m.b Xk]
        @debug("rect $rect")

        # gaussian elimination
        c, r = size(rect)

        rect[outgoing,:] /= rect[outgoing,end]

        for i in 1:c
        	if i != outgoing
        		rect[i,:] -= rect[outgoing,:] * rect[i,end]
        	end
        end

        @debug("rect $rect")

        m.b = rect[:,end-1]

    	invB[:,2:end] = rect[:,1:end-2]

    	m.basis[outgoing], m.nonbasis[incoming] = m.nonbasis[incoming], m.basis[outgoing]

        iter += 1
        if iter >= m.maxiter
            m.status = :UserLimit
            break
        end
    end

    m.objval = -m.b[1]

    m.solution = zeros(length(nonbasis_orig))

    @debug("basis $(m.basis)")
    @debug("nonbasis_orig $nonbasis_orig")
    @debug("b $(m.b)")
    @debug("objval $(m.objval)")
    @debug("solution $(m.solution)")

    for (i, v) in enumerate(nonbasis_orig)
        j = findfirst(x -> x == v, m.basis)
        if j > 0
            m.solution[i] = m.b[j]
        end
    end

    @debug("solution $(m.solution)")
end

function status(m::DrsMathProgModel)
    @debug("status")
    m.status
end

function getreducedcosts(m::DrsMathProgModel)
    @debug("getreducedcosts")
end

function getconstrduals(m::DrsMathProgModel)
    @debug("getconstrduals")
end

function getobjval(m::DrsMathProgModel)
    @debug("getobjval")
    m.objval
end

function getsolution(m::DrsMathProgModel)
    @debug("getsolution")
    m.solution
end

end

max_stars_repo_name: abogushov/lplib

max_issues_repo_name: abogushov/lplib

max_forks_repo_name: abogushov/lplib


text:
# Inception
function presentation_inception()
    f = conventional_inception()
    cache = AutoTM.Experiments.CPU_CACHE
    suffix = "study"

    # Custom ylabel
    ylabel = "Slowdown\\\\\\textbf{Lower is Better}"
    xlabel = "Dram Limit (GB)\\\\\\textbf{Lower is Better}"

    # Performance line graph
    only_dram = pgf_plot_performance(f, cache, suffix;
        file = joinpath(FIGDIR, "presentation_inception_dram.pdf"),
        formulations = ("synchronous",),
        filter_function = (x, y) -> ([first(x)], [first(y)]),
        ylabel = ylabel,
        xlabel = xlabel,
    )

    only_pmem = pgf_plot_performance(f, cache, suffix;
        file = joinpath(FIGDIR, "presentation_inception_pmm.pdf"),
        formulations = ("synchronous",),
        filter_function = (x, y) -> ([last(x)], [last(y)]),
        ylabel = ylabel,
        xlabel = xlabel,
    )

    all = pgf_plot_performance(f, cache, suffix;
        file = joinpath(FIGDIR, "presentation_inception_all.pdf"),
        formulations = ("synchronous",),
        ylabel = ylabel,
        xlabel = xlabel,
    )

    return (
        only_dram = only_dram,
        only_pmem = only_pmem,
        all = all,
    )
end

function presentation_2lm()
    fns = (large_vgg(), large_inception(), large_resnet(), large_densenet())
    cache = AutoTM.Experiments.CPU_CACHE
    cache_2lm = "nocache"
    formulations = ("static", "synchronous")

    pgf_large_performance(
        fns,
        cache,
        cache_2lm;
        formulations = formulations,
        ylabel = "Speedup over 2LM\\\\\\textbf{Higher is Better}",
        file = joinpath(FIGDIR, "presentation_large_performance.pdf")
    )
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
"""
    SOCtoRSOCBridge{T} <: Bridges.Variable.SetMapBridge{T,MOI.RotatedSecondOrderCone,MOI.SecondOrderCone}

Same transformation as [`MOI.Bridges.Constraint.SOCtoRSOCBridge`](@ref).
"""
struct SOCtoRSOCBridge{T} <:
       SetMapBridge{T,MOI.RotatedSecondOrderCone,MOI.SecondOrderCone}
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{
        MOI.VectorOfVariables,
        MOI.RotatedSecondOrderCone,
    }
end

"""
    RSOCtoSOCBridge{T} <: Bridges.Variable.SetMapBridge{T,MOI.SecondOrderCone,MOI.RotatedSecondOrderCone}

Same transformation as [`MOI.Bridges.Constraint.RSOCtoSOCBridge`](@ref).
"""
struct RSOCtoSOCBridge{T} <:
       SetMapBridge{T,MOI.SecondOrderCone,MOI.RotatedSecondOrderCone}
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.SecondOrderCone}
end


text:
# ---
# title: CVI-ICVI Comparison
# id: cvi_icvi
# cover: ../assets/logo.png
# date: 2021-12-6
# author: "[Sasha Petrenko](https://github.com/AP6YC)"
# julia: 1.6
# description: This demo illustrates how to use incremental training methods vs. batch training for all ART modules. This demo also demonstrates how the final results of each CVI and ICVI are equivalent.
# ---

# ## Overview

# This demo is a simple example of how to use CVIs incrementally and in batch to demonstrate that the end results are the same.
# Here, we load a simple dataset and run a basic clustering algorithm to prescribe a set of clusters to the features.
# In the incremental case, we will take advantage of the fact that we can compute a criterion value at every step by running the ICVI alongside an online clustering algorithm.
# This simple example demonstrates the usage of a single CVI/ICVI, but it may be substituted for any other CVI/ICVI in the `ClusterValidityIndices.jl` package.

# ## Data Setup

# First, we must load all of our dependencies.
# We will load the `ClusterValidityIndices.jl` along with some data utilities and the Julia `Clustering.jl` package to cluster that data.
using ClusterValidityIndices    # CVI/ICVI
using AdaptiveResonance         # DDVFA
using MLDatasets                # Iris dataset
using MLDataUtils               # Shuffling and splitting
using Printf                    # Formatted number printing
using Plots

# We will download the Iris dataset for its small size and benchmark use for clustering algorithms.
Iris.download(i_accept_the_terms_of_use=true)
features, labels = Iris.features(), Iris.labels()

# Because the MLDatasets package gives us Iris labels as strings, we will use the `MLDataUtils.convertlabel` method with the `MLLabelUtils.LabelEnc.Indices` type to get a list of integers representing each class:
labels = convertlabel(LabelEnc.Indices{Int}, labels)
unique(labels)

# ## CVI/ICVI Setup

# Because CVI/ICVIs only differ in their evaluation mode, we will use the same default constructor for both of our objects.
# Here, we will use the same type of CVI for both objects to verify that the results of both are the same at the final iteration.

## Create both CVI objects, using one incrementally as an ICVI
icvi = CH()
cvi = CH()

# ## Online Clustering

# ### Adaptive Resonance Theory Algorithms

# Adaptive Resonance Theory (ART) is a neurocognitive theory that is the basis of a class of online clustering algorithms.
# Because these clustering algorithms run online, we can both cluster and compute a new criterion value at every step.
# For more on these ART algorithms, see [`AdaptiveResonance.jl`](https://github.com/AP6YC/AdaptiveResonance.jl).

## Create a Distributed Dual-Vigilance Fuzzy ART (DDVFA) module with default options
art = DDVFA()
typeof(art)

# Because we are streaming clustering, we must setup the internal data setup of the DDVFA module.
# This is akin to doing some data preprocessing and communicating the dimension of the data, bounds, etc. to the module beforehand.

## Setup the data configuration for the module
data_setup!(art.config, features)
## Verify that the data is setup
art.config.setup

# ### Iteration

# We can now cluster and get the criterion values online
# We will do this by creating an ICVI object, setting up containers for the iterations, and then iterating.

## Setup the online/streaming clustering
n_samples = length(labels)          # Number of samples
c_labels = zeros(Int, n_samples)    # Clustering labels
criterion_values = zeros(n_samples) # ICVI outputs

## Iterate over all samples
for ix = 1:n_samples
    ## Extract one sample
    sample = features[:, ix]
    ## Cluster the sample online
    c_labels[ix] = train!(art, sample)
    ## Get the new criterion value (ICVI output)
    criterion_values[ix] = get_icvi!(icvi, sample, c_labels[ix])
end

## See the list of criterion values
criterion_values

# ## Batch Evaluation

# In batch mode, we will use the sample features and prescribed cluster labels as before to verify that the criterion values are the same at the last iteration.

## Get the final criterion value in batch
batch_criterion_value = get_cvi!(cvi, features, c_labels)

# ## Comparison

# Now we check that the two produce the same results in the end.

## Print the batch result and the final result of the incremental variant
@printf "Batch criterion value: %.4f\n" batch_criterion_value
@printf "Final incremental criterion value: %.4f\n" criterion_values[end]


text:
# SPDX-License-Identifier: MIT
# Author: Pierre Lairez

@genpredicate nogeneric function closest(p :: 3, q :: 3, a :: 3)
    qp = q - p
    pa = p - a
    qa = q - a

    Codegen.group!(qp...)
    Codegen.group!(pa..., qa...)

    inp(qp, pa+qa)
end
@doc """
    closest(p :: 3, q :: 3, a :: 3) -> Int

Return 1 if `a` is closer to `p` than to `q`.
Return –1 if `a` is closer to `q` than to `p`.
Return 0 is `a` is equaly close to both.
""" closest(::NTuple{3, Float64},::NTuple{3, Float64},::NTuple{3, Float64})

@genpredicate function orient(p :: 3, q :: 3, r :: 3, a :: 3)
    pa = p - a
    qa = q - a
    ra = r - a

    Codegen.group!(pa...)
    Codegen.group!(qa...)
    Codegen.group!(ra...)

    det(pa..., qa..., ra...)
end
@doc """
    orient(p :: 3, q :: 3, r :: 3, a :: 3) -> Int

Consider the oriented plane on which the triangle `pqr` is positively oriented.

* Return 1 if `a` is below this plane.
* Return –1 if `a` is above this plane.
* Return 0 if `a` lies on this plane.

""" orient(::NTuple{3, Float64},::NTuple{3, Float64},::NTuple{3, Float64},::NTuple{3, Float64})

@genpredicate function insphere(p :: 3, q :: 3, r :: 3, s :: 3, a :: 3)
    p = p - a
    q = q - a
    r = r - a
    s = s - a

    Codegen.group!(p..., q..., r..., s...)

    det(p..., abs2(p), q..., abs2(q), r..., abs2(r), s..., abs2(s))
end
@doc """
    insphere(p :: 3, q :: 3, r :: 3, s :: 3, a :: 3)

* Return 1 if `a` is inside the circumscribed sphere defined by the four points `p`, `q`, `r` and `s`.
* Return –1 if `a` is outside.
* Return 0 is `a` lies on the sphere or if the four points are coplanar.

""" insphere(::NTuple{3, Float64},::NTuple{3, Float64},::NTuple{3, Float64},::NTuple{3, Float64},::NTuple{3, Float64})


text:
using StatsPlots
using Query
using DataFrames
using Distributions
using LaTeXStrings
include("load_data.jl")

using Plots: PX_PER_INCH
const LINEWIDTH = 6.75*PX_PER_INCH
const COLWIDTH = 3.25*PX_PER_INCH
height(width, ratio=16/9) = width / ratio

dpl = 6.75 * 100
dph = dpl / (16/9)

plot_kwargs = (
    size=(1000, 400),
    markersize=8,
    layout=(2, 4),
    tickfontsize=8,
    guidefontsize=9,
    titlefontsize=10,
    legendfontsize=8,
    legendtitlefontsize=9,
)

####################################################################################################
# Performance: Alg, Order, Smoothing, Sigma
####################################################################################################
for probname in unique(df.probname)
    _df = df |>
        @filter((_.steprule .== :standard)) |>
        @filter((_.probname .== probname)) |>
        @filter((_.L2 .< 1e2)) |>
        DataFrame

    xlims, ylims = Dict(
        "logistic" => ((5e0, 1e6), (1e-13, 1e1),),
        "lotkavolterra" => ((1e2, 3e6), (1e-12, 1e2)),
        "fitzhughnagumo" => ((1e2, 3e6), (1e-13, 1e2)),
        "rigidbody" => ((4e1, 1.5e6), (5e-14, 1e2)),
        "brusselator" => ((4e1, 2.5e6), (1e-12, 2e2),),
    )[probname]

    legendplot = plot(
        repeat([missing], 1, 10),
        label=["EKF-IWP1" "EKF-IWP2" "EKF-IWP3" "EKF-IWP4" "EKF-IWP5" "EKS-IWP1" "EKS-IWP2" "EKS-IWP3" "EKS-IWP4" "EKS-IWP5"],
        color=[1 2 3 4 5 1 2 3 4 5],
        linestyle=hcat(repeat([:solid], 5)..., repeat([:dash], 5)...),
        marker=hcat(repeat([:o], 5)..., repeat([:star], 5)...),
        grid=false, xlims=(20,3), showaxis=false,
        legend=:left,
    );

    # plots = []
    i = 0
    for alg in ["EKF0", "EKF1"],
        sigma in [:fixedMLE, :dynamic, :fixedMLEMV,]
        i+=1

        if sigma in (:fixedMLEMV,) && alg == "EKF1"
            sigma = :dynamicMV
            alg = "EKF0"
        end

        sigmaname = Dict(
            :fixedMLE => "fixed",
            :dynamic => "TV",
            :fixedMLEMV => "fixed-MV",
            :dynamicMV => "TV-MV",
        )[sigma]

        p = _df |>
            @filter((_.sigmarule .== sigma)) |>
            @filter((_.alg .== alg)) |>
            @df StatsPlots.plot(
                :num_evals, :L2,
                group=(:q_str, :smooth_str),
                color=:q,
                linestyle=[:solid, :dash, :dot, :dashdot][:smooth.+1],
                markershape=[:o, :utriangle][:smooth.+1],
                xscale=:log10, yscale=:log10,
                xlabel = (alg == "EKF1" || sigma in (:dynamicMV,)) ? L"\#Evaluations" : " ",
                ylabel = sigma == :fixedMLE ? "L2 Error" : "",
                title="EKF$(alg[end])/EKS$(alg[end]), Γ=$sigmaname",
                legend=false,
            );
        plot!(p, xlims=xlims, ylims=ylims,)
        savefig(p, "/home/nath/PhD-Projects/ode-filter-uq-paper/figures/$(probname)_performance_$i.tex")

    end
    l = @layout [
        [a b c
         d e f] g{0.13w}
    ]
    p = plot(plots..., legendplot;
         xlims=xlims, ylims=ylims,
         plot_kwargs...,
         layout=l,
         );
    pgfsave("figures/$(probname)_performance.pdf")
    savefig("figures/$(probname)_performance.pdf")
end

####################################################################################################
# Calibration
####################################################################################################
for probname in unique(df.probname)
    _df = df |>
        @filter((_.steprule .== :standard)) |>
        @filter((_.probname .== probname)) |>
        @filter((_.L2 .< 1e2)) |>
        DataFrame

    xlims, ylims = Dict(
        "logistic" => ((1e-13, 1e1), (1e-7, 1e3),),
        "lotkavolterra" => ((1e-12, 1e2), (5e-9, 1.1e3)),
        "fitzhughnagumo" => ((1e-13, 1e2), (1e-8, 1e4)),
        "rigidbody" => ((5e-14, 1e2), (8e-7, 2e4)),
        "brusselator" => ((1e-12, 2e2), (5e-9, 1e3),),
    )[probname]

    legendplot = plot(
        repeat([missing], 1, 12),
        label=["EKF-IWP1" "EKF-IWP2" "EKF-IWP3" "EKF-IWP4" "EKF-IWP5" "EKS-IWP1" "EKS-IWP2" "EKS-IWP3" "EKS-IWP4" "EKS-IWP5" "E[χ²(2)]" "99% interval"],
        color=[1 2 3 4 5 1 2 3 4 5 :black :black],
        linestyle=hcat(repeat([:solid], 5)..., repeat([:dash], 5)..., :dash, :dot),
        marker=hcat(repeat([:o], 5)..., repeat([:star], 5)..., :none, :none),
        grid=false, xlims=(20,3), showaxis=false,
        legend=:left,
    )

    plots = []
    for alg in ["EKF0", "EKF1"],
        sigma in [:fixedMLE, :dynamic, :fixedMLEMV]

        if sigma in (:fixedMLEMV,) && alg == "EKF1"
            sigma = :dynamicMV
            alg = "EKF0"
        end

        sigmaname = Dict(
            :fixedMLE => "fixed",
            :dynamic => "TV",
            :fixedMLEMV => "fixed-MV",
            :dynamicMV => "TV-MV",
        )[sigma]

        p = plot()
        _p = 1-0.99
        dist = Chisq(2)
        m, q_low, q_high = mean(dist), quantile(dist, 1-(_p/2)), quantile(dist, _p/2)
        hline!(p, [m], color=:black, linestyle=:dash, label="")
        hline!(p, [q_low, q_high], color=:black, linestyle=:dot, label="")

        _df |>
            @filter((_.sigmarule .== sigma)) |>
            @filter((_.alg .== alg)) |>
            @df StatsPlots.plot!(p,
                # :num_evals, :Χ²,
                :L2, :Χ²,
                group=(:q_str, :smooth_str),
                color=:q,
                linestyle=[:solid, :dash, :dot, :dashdot][:smooth.+1],
                markershape=[:o, :star][:smooth.+1],
                xscale=:log10, yscale=:log10,
                xlabel = (alg == "EKF1" || sigma == :dynamicMV) ? "L2 Error" : " ",
                ylabel = sigma == :fixedMLE ? "χ²-statistic" : "",
                title="EKF$(alg[end])/EKS$(alg[end]), Γ=$sigmaname",
                legend=false,
            )
        push!(plots, p)
    end
    l = @layout [
    [a b c
     d e f] g{0.13w}
    ]
    plot(plots..., legendplot;
         xlims=xlims, ylims=ylims,
         plot_kwargs...,
         layout=l,
         )
    savefig("figures/$(probname)_calibration.pdf")
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
# data type and constructors to mimic some internal definitions of julia

mutable struct NewModule
    name::Symbol
    parent::Union{NewModule,Void}
    constants::Dict{Any,Any}
end

abstract type NewType end

struct NewBottomType <: NewType end

struct NewTypeName
    name::Symbol
    modul::NewModule
end

NewTypeBounds = Union{NewType,Any}
mutable struct NewTypeVar
    name::Symbol
    lb::NewTypeBounds
    ub::NewTypeBounds
end
NewTypeBounds = Union{NewType,NewTypeVar}

mutable struct NewDataType <: NewType
    name::NewTypeName
    super::Any # NewType
    parameters::Vector
    abstr::Bool
    mutable::Bool
    isleaftype::Bool
    function NewDataType(name, super, parameters, abstr, mutable, isleaftype)
        dt = new(name, super, parameters, abstr, mutable, isleaftype)
        register!(dt)
    end
end

struct DataTypeKey
    name::Symbol
    parameters::Vector
end

struct NewUnionAll <: NewType
    var::NewTypeVar
    body::NewType
end

struct NewUnion <: NewType
    a
    b
    NewUnion(a::NewType, b::NewType, ::Bool) = new(a, b)
end

mutable struct NewMethod
    name::Symbol
    modul::NewModule
    file::Symbol
    line::Int32
    sig::Union{NewDataType,NewUnionAll}
    nargs::Int32
    pure::Bool
end

#mutable struct NewMethodInstance
#end

mutable struct NewTypeMapEntry
    next::Union{NewTypeMapEntry, Void}
    sig::Tuple
    func::NewMethod
    isleafsig::Bool
    issimplesig::Bool
    va::Bool
end

mutable struct NewMethodTable
    name::Symbol
    defs::Union{NewTypeMapEntry,Void}
    max_args::Int64
    modul::NewModule
end

mutable struct NewMethodList
    ms::Array{NewMethod,1}
    mt::NewMethodTable
    NewMethodList(ms::Array{NewMethod,1}, mt::NewMethodTable) = new(ms, mt)
end

struct FreeVariables
    d::Dict{NewTypeVar,Any} # Array{Union{Void,NewTypeVar,NewType,Bitstype}}
end


text:
module DivConformingTestsSeq
using PartitionedArrays
include("../DivConformingTests.jl")
prun(DivConformingTests.main,sequential,2)
end # module


text:
function file_transfer_test()
    @testset "File transfer" begin
        srcDir = "/Users/lutz/Documents/temp";
        tgDir = "/nas/longleaf/home/lhendri/Documents/temp"

        test_header("Common base dir")
        # Relative path
        srcPath = "abc/def.txt";
        baseDirV = ["/base/one", "/two/base", "/base/three"];
        outDirV = common_base_dir(srcPath, baseDirV);
        for j = 1 : length(baseDirV)
            @test isequal(outDirV[j],  joinpath(baseDirV[j], srcPath))
        end
        # Absolute path
        srcPathAbs = joinpath(baseDirV[2], srcPath);
        outDir2V = common_base_dir(srcPathAbs, baseDirV);
        @test all(isequal.(outDirV, outDir2V))

        # Remote and local path using Computers as inputs
        localDir, remoteDir = remote_and_local_path(srcPath,  tgComp = :longleaf);
        @test localDir == joinpath(home_dir(FilesLH.defaultLocal), srcPath)
        @test remoteDir == joinpath(home_dir(:longleaf), srcPath)

        # The same with full path
        localDir2, remoteDir2 = remote_and_local_path(
            joinpath(homedir(), srcPath),  tgComp = :longleaf);
        @test isequal(localDir2, localDir)
        @test isequal(remoteDir, remoteDir2)

        test_header("Remote copy");
        localPath = joinpath(srcDir, "temp1.txt");
        remotePath = joinpath(tgDir, "temp1.txt");
        sshStr = "lhendri@longleaf.unc.edu";
        cmdUp = remote_copy(localPath, remotePath, sshStr, :up, trialRun = true);
        @test isa(cmdUp, Cmd)
        cmdDown = remote_copy(remotePath, localPath, sshStr, :up, trialRun = true);
        @test isa(cmdDown, Cmd)

        # This creates unreasonable amounts of output. Redirect or change dir +++
        test_header("Rsync")
        rsCmd = rsync_command(srcDir, tgDir; trialRun = true, verbose = false,
            excludeV = "abc/test.txt");
        @test isa(rsCmd, Cmd);
        rsCmd = rsync_command(srcDir, tgDir; trialRun = true, verbose = false,
            excludeV = ["abc/test.txt", "def/ghi"]);
        @test isa(rsCmd, Cmd);

        rsCmd = FilesLH.rsync_dir(srcDir, tgDir; trialRun = true, verbose = false);
        @test isa(rsCmd, Cmd);

        test_header("Git")
        repoDir, _ = splitdir(@__DIR__);
        @test git_upload_dir(repoDir, trialRun = true);
    end
end

# --------------------


text:
using Test, SeisRequests
using Dates: Second, Millisecond

@testset "Utils" begin
    @testset "seconds_milliseconds" begin
        @test SeisRequests.seconds_milliseconds(5.9019) == Second(5) + Millisecond(901)
    end

    @testset "split_channel_code" begin
        @test SeisRequests.split_channel_code("AB,CD.*..BH?") ==
            (network="AB,CD", station="*", location="", channel="BH?")
        @test_throws ArgumentError SeisRequests.split_channel_code("A.B.C.D.")
        @test_throws ArgumentError SeisRequests.split_channel_code("A.B.C")
    end
end


text:
export STDP

abstract type LearningRule end

"""
# STDP learning rule
Spike Timing Dependent Plasticity learning rule. See: http://www.scholarpedia.org/article/Spike-timing_dependent_plasticity
# Arguments
...
- `trace_aditive::Bool`: if true performs a "all-to-all interaction" and else performs a "nearest-neighbor interaction".
"""
@with_kw mutable struct STDP <: LearningRule
    A₊::Real
    A₋::Real
    τ₊::Real
    τ₋::Real
    min_weight::Union{Real,AbstractMatrix} = -Inf
    max_weight::Union{Real,AbstractMatrix} = Inf
    traces_additive::Bool = false
    update_rule::UpdateRule = RegularUpdate()
    e₊::Union{AbstractArray,Nothing} = nothing
    e₋::Union{AbstractArray,Nothing} = nothing
end

function add_group!(learning_rule::STDP, N::Int, batch_size::Int)
    if isnothing(learning_rule.e₊)
        learning_rule.e₊ = zeros(Float32, batch_size, 0)
        learning_rule.e₋ = zeros(Float32, batch_size, 0)
    end
    learning_rule.e₊ = pad1D(learning_rule.e₊, N)
    learning_rule.e₋ = pad1D(learning_rule.e₋, N)
end

function reset!(learning_rule::STDP)
    fill!(learning_rule.e₊, 0)
    fill!(learning_rule.e₋, 0)
end


text:
import Dates
using Base: Semaphore
using Base.Threads: Atomic, @spawn, Condition
using Dates: Period

function catch_stack_to_string(catch_stack)
    buf = IOBuffer()
    Base.display_error(buf, catch_stack)
    return String(take!(buf))
end
function restore_callsite_source_position!(expr, src)
    # because the logging macros call functions in this package, the source file+line are
    # incorrectly attributed to the wrong location, so we explicitly override the source
    # with the original value
    expr.args[1].args[2] = src
    return expr
end
macro error_with_catch_stack(msg, exs...)
    return restore_callsite_source_position!(esc(:(
        $Base.@error string($msg, "\n", catch_stack_to_string($Base.catch_stack())) $(exs...)
        )), __source__)
end

"""
    @spawn_with_error_log expr
    @spawn_with_error_log "..error msg.." expr
Exactly like `@spawn`, except that it wraps `expr` in a try/catch block that will print any
exceptions that are thrown from the `expr` to stderr, via `@error`. You can
optionally provide an error message that will be printed before the exception is displayed.
This is useful if you need to spawn a "background task" whose result will never be
waited-on nor fetched-from.
"""
macro spawn_with_error_log(expr)
    _spawn_with_error_log_expr(expr)
end
macro spawn_with_error_log(message, expr)
    _spawn_with_error_log_expr(expr, message)
end
function _spawn_with_error_log_expr(expr, message = "")
    e = gensym("err")
    return esc(
        quote
            $Base.Threads.@spawn try
                $(expr)
            catch $e
                $TransactionLogging.@error_with_catch_stack "@spawn_with_error_log failed:"
                rethrow()
            end
        end
    )
end

"""
    PeriodicTask
This structure is a wrapper around background periodic Task and can be used to inspect the
state of the task itself and to safely terminate the background periodic task by signalling
via `should_terminate`.
"""
struct PeriodicTask
    # Name of the periodic task. Attached to error logs for debuggability.
    name::String

    # Specifies how often the underlying periodic task should be run.
    period::Period

    # When set to true, the underlying periodic task will terminate before next
    # iteration.
    should_terminate::Atomic{Bool}

    # This is used by the task to wait until `period` elapses and is used
    # in stop_periodic_task! to quickly wake-up and terminate sleeping task.
    condition::Condition

    # The underlying periodic task itself.
    task::Task
end

"""
    @spawn_periodic_task period expr name
Run `expr` once every `period` and returns `PeriodicTask` that will carry out this logic. The task
can be terminated by calling `stop_periodic_task!`. Optional `name` can be specified and
this will be attached to error logs and (eventually) metrics for easier debuggability.
`period` must be a `Dates.Period`.
# Examples
```julia
import Dates
disk_stuff = DiskStuff()
my_task = @spawn_periodic_task Dates.Seconds(30) dump_some_stuff_to_disk(disk_stuff) "DiskDumper"
# ... do some stuff ...
stop_periodic_task!(my_task)
istaskfailed(my_task) && throw(SystemError("Periodic task has failed!!"))
```
"""
macro spawn_periodic_task(period, expr, name="Unnamed")
    # TODO(janrous): add number of iteratons, number of failures and last successful
    # iteration timestamp once metrics with labels are available.
    return quote
        should_terminate = Atomic{Bool}(false)
        cond = Condition()
        # TODO(janrous): we can improve the timing precision by calculating how
        # much time has elapsed since last_execution when we are about to enter
        # sleep and subtract that from period.
        task = @spawn begin
            while !should_terminate[]
               Base.@lock cond begin
                   # Spawn notifier to wake us up when it's time to run.
                    @spawn begin
                        sleep($(esc(period)))
                        Base.@lock cond notify(cond)
                    end
                    wait(cond)
                    should_terminate[] && continue
                    try
                        $(esc(expr))
                    catch err
                        @error_with_catch_stack "$($(esc(name))): periodic task failed"
                    end
                end
            end
        end
        # TODO: if name is not given, use module:lineno of the caller
        PeriodicTask($(esc(name)), $(esc(period)), should_terminate, cond, task)
    end
end

macro spawn_sticky_periodic_task(name, period, expr)
    return quote
        should_terminate = Atomic{Bool}(false)
        cond = Condition()
        # This version uses @async instead of @spawn to ensure that tasks are run
        # on the current thread in a "sticky" mode. This is intended to be used
        # with the high priority julia patch that moves all generic task handling
        # from the thread that calls this macro.
        task = @async begin
            @info "Scheduled sticky periodic task $($(esc(name))) on thread $(Threads.threadid())"
            while !should_terminate[]
               Base.@lock cond begin
                   # Spawn notifier to wake us up when it's time to run.
                    @async begin
                        sleep($(esc(period)))
                        Base.@lock cond notify(cond)
                    end
                    wait(cond)
                    should_terminate[] && continue
                    try
                        $(esc(expr))
                    catch err
                        @error_with_catch_stack "$($(esc(name))): periodic task failed"
                    end
                end
            end
        end
        PeriodicTask($(esc(name)), $(esc(period)), should_terminate, cond, task)
    end
end

"""
    stop_periodic_task!(task::PeriodicTask)
Triggers termination of the periodic task.
"""
function stop_periodic_task!(task::PeriodicTask)
    task.should_terminate[] = true
    Base.@lock task.condition notify(task.condition)
    wait(task.task)
    return task
end

# Reflection methods for the inner ::Task struct.
Base.istaskdone(t::PeriodicTask) = istaskdone(t.task)
Base.istaskfailed(t::PeriodicTask) = istaskfailed(t.task)
Base.istaskstarted(t::PeriodicTask) = istaskstarted(t.task)


text:
using CausalityTools,StatsBase, DelayEmbeddings, LinearAlgebra, Statistics
import DelayEmbeddings.Neighborhood: Theiler, bulkisearch, NeighborNumber

export crossmap

"""
    crossmap(x, y, d, τ; r = 0, correspondence_measure = Statistics.cor) → Float64
    crossmap(x, y, d, τ, bootstrap_method::Symbol; r = 0, correspondence_measure = Statistics.cor,
        method = :segment, L = ceil(Int, (length(x)-d*τ)*0.2), nreps = 100) → Vector{Float64}

Compute the cross mapping [^Sugihara2012] between `x` and `y`, which is the correspondence (computed using
`correspondence measure`) between the values ``y(t)`` and the cross-map estimated values ``ỹ(t) | M_x``.
Returns the correspondence between original and cross mapped values (the default is
`ρ = correspondence_measure(y(t), ỹ(t) | M_x)`).

Here, ``y(t)`` are the raw values of the time series `y`, and ``ỹ(t)`` are the predicted values
computed from the out-of-sample embedding ``M_X`` constructed from the time series `x` with
embedding dimension `d` and embedding lag `τ`.

The Theiler window `r` indicates how many temporal neighbors of the predictee is to be excluded
during the nearest neighbors search (the default `r = 0` excludes only the predictee itself, while
`r = 2` excludes the point itself plus its two nearest neighbors in time).

If `bootstrap_method` is specified, then `nreps` different bootstrapped estimates of
`correspondence_measure(y(t), ỹ(t) | M_x)` are returned. The following bootstrap methods are available:

- `bootstrap_method = :random` selects training sets of length `L` consisting of randomly selected
    points from the embedding ``M_x``  (time ordering does not matter). This is method 3 from Luo
    et al. (2015)[^Luo2015], which critiqued the original Sugihara et al. methodology.
- `bootstrap_method = :segment` selects training sets consisting of time-contiguous segments
    (each of lenght `L`) of embedding vectors in ``M_x`` (time ordering matters). This is
    method 2 from Luo et al. (2015)[^Luo2015].

[^Sugihara2012]: Sugihara, George, et al. "Detecting causality in complex ecosystems." Science (2012): 1227079.[http://science.sciencemag.org/content/early/2012/09/19/science.1227079](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)
[^Luo2015]: "Questionable causality: Cosmic rays to temperature." Proceedings of the National Academy of Sciences Aug 2015, 112 (34) E4638-E4639; DOI: 10.1073/pnas.1510571112 Ming Luo, Holger Kantz, Ngar-Cheung Lau, Wenwen Huang, Yu Zhou
"""
function crossmap(x, y, d, τ; correspondence_measure = Statistics.cor, r = 0)
    # Embed `x`
    Mₓ = crossmapembed(x, d, τ, CCMEmbedding())

    # Compute `ỹ | Mₓ` and return `correspondence_measure(ỹ | Mₓ, y).
    return crossmap_basic(Mₓ, y, d, τ, correspondence_measure = correspondence_measure, r = r)
end

function crossmap(x, y, d, τ, bootstrap_method::Symbol;
        L = ceil(Int, (length(x) - d * τ) * 0.2), nreps = 100,
        r = 0, correspondence_measure = Statistics.cor)

    # Embed `x`
    @assert length(x) == length(y)
    @assert length(x) - d * τ > 0
    Mₓ = crossmapembed(x, d, τ, CCMEmbedding());

    # Compute `ỹ | Mₓ` and return `nreps` values for `correspondence_measure(ỹ | Mₓ, y),
    # each value being computed on a random library (sequential or not) consisting of `L` points.
    return crossmap_bootstrap(Mₓ, y, d, τ, bootstrap_method; L = L, nreps = nreps, r = r,
        correspondence_measure = correspondence_measure)
end


text:
using VLMinorityGameKit
using BSON
using Plots
using DataFrames

# setup -
# path_to_game_world = "/Users/jeffreyvarner/Desktop/julia_work/VLMinorityGameKit.jl/test/data/GW-M10-S5-NA1001.bson"
path_to_game_world = "/Users/jeffreyvarner/Desktop/julia_work/VLMinorityGameKit.jl/test/data/GW-M10-S5-NA1001.bson"
path_to_results_dir = "/Users/jeffreyvarner/Desktop/julia_work/VLMinorityGameKit.jl/test/data/results/SIM-M10-S5-NA1001-L3.bson"
d = BSON.load(path_to_game_world)
gameWorld = d[:gameworld]

println("Building the game world - finished. Starting the simulation ...")

# setup the number of time steps -
number_of_samples = 100
numberOfTimeSteps = 140
numberOfTraders = 1001
price_array = Array{Float64,2}(undef, (numberOfTimeSteps), number_of_samples)
simulation_result_array = Array{DataFrame,1}(undef, number_of_samples)

for sample_index = 1:number_of_samples
    result = execute(gameWorld, numberOfTimeSteps; liquidity=3.0 * numberOfTraders);
    simulation_result_array[sample_index] = result.market_table
    println("Completed sample = $(sample_index) ...")
end

# package -
simulation_dictionary = Dict()
simulation_dictionary["number_of_samples"] = number_of_samples
simulation_dictionary["number_of_timesteps"] = numberOfTimeSteps
simulation_dictionary["number_of_traders"] = numberOfTraders
simulation_dictionary["memory"] = 10
simulation_dictionary["strategies_per_agent"] = 5
simulation_dictionary["liquidity"] = 3.0 * numberOfTraders
simulation_dictionary["simulation_result_array"] = simulation_result_array

# dump -
bson(path_to_results_dir, data=simulation_dictionary)

# result = simulate(gameWorld, numberOfTimeSteps; liquidity=10.0 * numberOfTraders);
println("Simulation is finished .... enjoy.")


text:
# TODO:  TO FIX !!!

using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using SparseArrays, LinearAlgebra
using ViewerGL
GL = ViewerGL

function randomshapes()
	V,EV,FV = Lar.randomcubes(8, .8)
	model3d = V,EV,FV
	VV = [[k] for k=1:size(V,2)]

	Sigma =  Lar.spaceindex((V,FV))
	GL.VIEW( GL.numbering(.15)((V,[VV, EV, FV])) )
	V,EVs = Lar.biconnectedComponent((V,EV::Lar.Cells)) # 2-connected components (H & T)
	GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,12,0.99));

	W,bases,coboundaries = Lar.chaincomplex(V,FV,EV)
	EV, FV, CV = bases
	copEV, copFE, copCF = coboundaries # OK

#
EV = Lar.cop2lar(copEV)
FE = [findnz(copFE[k,:])[1] for k=1:size(copFE,1)]
FV = [union([EV[e] for e in FE[f]]...) for f=1:length(FE)]
W = convert(Lar.Points, W)
WW = [[k] for k=1:size(W,2)]

GL.VIEW(GL.numbering(.2)((W,[WW,EV]) ));
#

	W = convert(Lar.Points, W')
	triangulated_faces = Lar.triangulate(W, [copEV, copFE])
	FVs = convert(Array{Lar.Cells}, Lar.triangulated_faces)

	Plasm.viewcolor(V::Lar.Points, FVs::Array{Lar.Cells})
end

randomshapes()