text:
###############################################################################
#                                                                             #
#  This file is part of the julia module for Multi Objective Optimization     #
#  (c) Copyright 2017 by Aritra Pal                                           #
#                                                                             #
#  Permission is granted for academic research use.  For other uses,          #
#  contact the authors for licensing options.                                 #
#                                                                             #
#  Use at your own risk. I make no guarantees about the correctness or        #
#  usefulness of this code.                                                   #
#                                                                             #
###############################################################################

import Modof, Modolib, Modoalgos, Match, DataFrames, DataFramesMeta
@everywhere using Modof, Modolib, Modoalgos, Match, DataFrames, DataFramesMeta

params = Dict()
params[:algorithm] = "BBM"
params[:type_of_operation] = "lex_min"
params[:mip_emphasis] = true
params[:type_of_boost] = "NB"
params[:turn_off_cplex_heur] = true
params[:timelimit] = Inf
params[:total_threads] = 1

@everywhere function solve_bopinstance_exactly(params)
	if params[:read_instance] != :read_boscp_xavier
		instance, true_non_dom_sols = eval(params[:read_instance])(params[:read_instance_input])
		println("Solving Instance = $(params[:read_instance_input])")
	else
		instance, true_non_dom_sols = eval(params[:read_instance])(params[:read_instance_input]...)
		println("Solving Instance = $(params[:read_instance_input]...)")
	end

	if params[:total_threads] == 1
		println("Solving Instance on a Single Thread")
		master0::Float64 = time()
		non_dom_sols, rects, ip_solved = compute_exact_frontier(instance, params)
		master1::Float64 = time()
		println("Finished")
		non_dom_sols = select_non_dom_sols(non_dom_sols)
		slave0::Float64, slave1::Float64 = 0.0, 0.0
	else
		println("Decomposing Instance into $(params[:total_threads]) Subproblems")
		master0 = time()
		non_dom_sols_, rects_, ip_solved_ = Modoalgos.decompose_a_problem_into_smaller_problem(instance, params)
		println("Finished Decomposition")
		master1 = time()

		algorithm = " "
		if params[:algorithm] == "AUTO"
			algorithm = "AUTO"
			density::Float64 = length(findn(vec(instance.A)))/(size(instance.A)[1]*size(instance.A)[2])
			if density >= 0.42
				params[:algorithm] = "EBBM"
				params[:type_of_boost] = "O_SH"
				params[:turn_off_cplex_heur] = true
				params[:type_of_operation] = "lex_min"
				params[:mip_emphasis] = true
			else
				params[:algorithm] = "UDECM"
				params[:type_of_boost] = "NB"
				params[:turn_off_cplex_heur] = false
				params[:type_of_operation] = "dynamic"
				params[:mip_emphasis] = false
				if floor(abs(log10(size(instance.A)[1]/size(instance.A)[2]))) > 1.0
					params[:algorithm] = "EBBM"
				end
			end
		end

		println("Started Solving Subproblem")
		slave0 = time()
		non_dom_sols, rects, ip_solved = @match params[:algorithm] begin
			"UDECM" || "BDECM" => Modoalgos.ϵ_constraint_method(rects_[params[:thread_id]], instance, non_dom_sols_[params[:thread_id]], params)
			_ => Modoalgos.balanced_box_method(rects_[params[:thread_id]], instance, non_dom_sols_[params[:thread_id]], params)
		end
		slave1 = time()
		ip_solved += ip_solved_
		println("Finished Solving Subproblem")
		non_dom_sols = select_non_dom_sols(non_dom_sols)
		if algorithm == "AUTO"
			params[:algorithm] = "AUTO"
			params[:type_of_boost] = "NB"
			params[:turn_off_cplex_heur] = false
			params[:type_of_operation] = "lex_min"
			params[:mip_emphasis] = false
		end
	end

	filename::String = "../results/Experimental_Results_$(myid()).csv"
	cor_::Float64 = cor(instance.c1, instance.c2)
	sparsity::Float64 = 1-(length(findn(vec(instance.A)))/(size(instance.A)[1]*size(instance.A)[2]))

	if isfile(filename) == false
		results::DataFrame = DataFrame(Problem_Type=params[:problem_type], Instance_Type=params[:instance_type], Instance_ID=params[:instance_id], Objectives=2, Correlation=cor_, Variables=size(instance.A)[2], Constraints=size(instance.A)[1], Sparsity=sparsity, Algorithm=params[:algorithm], Type_of_Operation=params[:type_of_operation], MIP_Emphasis=params[:mip_emphasis], Type_of_Boost=params[:type_of_boost], Cplex_Heur=!params[:turn_off_cplex_heur], Parallelism=params[:parallelism], Total_Threads=params[:total_threads], Thread_ID=params[:thread_id], Master_Time=master1-master0, Slave_Time=slave1-slave0, IP_Solved=ip_solved, Num_Non_Dom_Pts=size(true_non_dom_sols)[1], Num_Non_Dom_Pts_Found=length(non_dom_sols))
	else
		results = readtable(filename)
		results[:Slave_Time] = float(results[:Slave_Time])
		push!(results, [params[:problem_type], params[:instance_type], params[:instance_id], 2, cor_, size(instance.A)[2], size(instance.A)[1], sparsity, params[:algorithm], params[:type_of_operation], params[:mip_emphasis], params[:type_of_boost], !params[:turn_off_cplex_heur], params[:parallelism], params[:total_threads], params[:thread_id], master1-master0, slave1-slave0, ip_solved, size(true_non_dom_sols)[1], length(non_dom_sols)])
	end
	writetable(filename, results)
end

@everywhere instance, true_non_dom_sols = read_bokp_xavier1("2KP50-1A")

###############################################################################
# Testing all Methods                                                         #
###############################################################################

@everywhere function count_true_non_dom_sols(true_::Array{Float64, 2}, non_dom_sols::Vector{BOPSolution})
	count = 0
	inds = Int64[]
	for i in 1:size(true_)[1]
		for j in 1:length(non_dom_sols)
			if true_[i, :] == [non_dom_sols[j].obj_val1, non_dom_sols[j].obj_val2]
				push!(inds, j)
				count += 1
				break
			end
		end
	end
	count, non_dom_sols[inds]
end

@everywhere function testing_all_methods(instance::BOBPInstance, true_non_dom_sols::Array{Float64, 2})
	params = Dict()
	params[:algorithm] = "UDECM"
	params[:type_of_operation] = "lex_min"
	params[:mip_emphasis] = false
	params[:type_of_boost] = "NB"
	params[:turn_off_cplex_heur] = false
	params[:timelimit] = Inf
	params[:total_threads] = 1

	println("Testing 1 Thread UDECM_LM_NB_HON")
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread BBM_LM_NB_HON")
	params[:algorithm] = "BBM"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread EBBM_LM_NB_HON")
	params[:algorithm] = "EBBM"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread UDECM_AO_NB_HON")
	params[:algorithm] = "UDECM"
	params[:type_of_operation] = "augmented_operation"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread EBBM_AO_NB_HON")
	params[:algorithm] = "EBBM"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread UDECM_D_NB_HON")
	params[:algorithm] = "UDECM"
	params[:type_of_operation] = "dynamic"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread EBBM_D_NB_HON")
	params[:algorithm] = "EBBM"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread EBBM_LM2_NB_HON")
	params[:algorithm] = "EBBM"
	params[:type_of_operation] = "lex_min"
	params[:mip_emphasis] = true
	params[:type_of_boost] = "NB"
	params[:turn_off_cplex_heur] = false
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread UDECM_D_B_HON")
	params[:algorithm] = "UDECM"
	params[:type_of_operation] = "dynamic"
	params[:mip_emphasis] = false
	params[:type_of_boost] = "B"
	params[:turn_off_cplex_heur] = true
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread UDECM_D_O_SH_HON")
	params[:type_of_boost] = "O_SH"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread EBBM_D_B_HON")
	params[:algorithm] = "EBBM"
	params[:type_of_boost] = "B"
	params[:turn_off_cplex_heur] = true
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread EBBM_D_O_SH_HON")
	params[:type_of_boost] = "O_SH"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end

	println("Testing 1 Thread AUTO")
	params[:algorithm] = "AUTO"
	non_dom_sols, rects = compute_exact_frontier(instance, params)
	non_dom_sols = select_non_dom_sols(non_dom_sols)
	count, non_dom_sols = count_true_non_dom_sols(true_non_dom_sols, non_dom_sols)
	if length(non_dom_sols) == size(true_non_dom_sols)[1]
		println("Test Successful")
	else
		println("Test Not Successful - Difference = $(size(true_non_dom_sols)[1]-count)")
	end
end

###############################################################################
# Experiments for All Methods                                                 #
###############################################################################

###############################################################################
## Single Thread                                                             ##
###############################################################################

@everywhere function generate_queue_for_detailed_experiments_on_all_instances()
	queue_ = Dict{Any, Any}[]

	solver_params = Dict()
	solver_params[:algorithm] = ["UDECM", "BBM", "EBBM", "UDECM", "UDECM", "EBBM", "EBBM", "EBBM", "UDECM", "UDECM", "EBBM", "EBBM", "EBBM", "EBBM", "EBBM", "EBBM", "AUTO", "AUTO", "AUTO"]
	solver_params[:type_of_operation] = ["lex_min", "lex_min", "lex_min", "augmented_operation", "dynamic", "augmented_operation", "dynamic", "lex_min", "dynamic", "dynamic", "dynamic", "dynamic", "lex_min", "lex_min", "lex_min", "lex_min", "lex_min", "lex_min", "lex_min"]
	solver_params[:mip_emphasis] = [false, false, false, false, false, false, false, true, false, false, false, false, true, true, true, true, false, false, false]
	solver_params[:type_of_boost] = ["NB", "NB", "NB", "NB", "NB", "NB", "NB", "NB", "B", "O_SH", "B", "O_SH", "B", "O_SH", "B", "O_SH", "NB", "NB", "NB"]
	solver_params[:turn_off_cplex_heur] = [false, false, false, false, false, false, false, false, false, false, false, false, false, false, true, true, false, false, false]
	solver_params[:total_threads] = [ones(16)..., 2, 3, 4]

	for i in [[100, 10], [200, 40], [600, 60]], j in ["a", "b", "c", "d"], k in 1:length(solver_params[:algorithm])
		for t in 1:solver_params[:total_threads][k]
			params = Dict()
			params[:read_instance] = :read_boscp_xavier
			params[:read_instance_input] = [i..., j]
			params[:problem_type] = "BOSCP"
			params[:instance_type] = "Xavier"
			params[:instance_id] = "SCP$(i[2])x$(i[1])_$(j)"
			params[:algorithm] = solver_params[:algorithm][k]
			params[:type_of_operation] = solver_params[:type_of_operation][k]
			params[:mip_emphasis] = solver_params[:mip_emphasis][k]
			params[:type_of_boost] = solver_params[:type_of_boost][k]
			params[:turn_off_cplex_heur] = solver_params[:turn_off_cplex_heur][k]
			params[:timelimit] = Inf
			params[:total_threads] = solver_params[:total_threads][k]
			params[:thread_id] = t
			push!(queue_, params)
		end
	end

	instances = readdir(Pkg.dir("Modolib/instances/BOKP/Xavier_Instances/non dominated points/"))
	ind = Int64[]
	for i in 1:length(instances)
		instances[i]= split(instances[i],".")[1]
		if contains(instances[i], "2KP")
			push!(ind, i)
		end
    end
    for i in 1:length(instances), k in 1:length(solver_params[:algorithm])
    	for t in 1:solver_params[:total_threads][k]
			params = Dict()
    		if i in ind
				params[:read_instance] = :read_bokp_xavier1
			else
				params[:read_instance] = :read_bokp_xavier2
			end
			params[:read_instance_input] = instances[i]
			params[:problem_type] = "BOKP"
			params[:instance_type] = "Xavier"
			params[:instance_id] = instances[i]
			params[:algorithm] = solver_params[:algorithm][k]
			params[:type_of_operation] = solver_params[:type_of_operation][k]
			params[:mip_emphasis] = solver_params[:mip_emphasis][k]
			params[:type_of_boost] = solver_params[:type_of_boost][k]
			params[:turn_off_cplex_heur] = solver_params[:turn_off_cplex_heur][k]
			params[:timelimit] = Inf
			params[:total_threads] = solver_params[:total_threads][k]
			params[:thread_id] = t
			push!(queue_, params)
		end
	end

	instances = readdir(Pkg.dir("Modolib/instances/BOSPP/Xavier_Instances/non dominated points/"))
	for i in 1:length(instances)
		instances[i]= split(instances[i],".")[1]
    end
    sort!(instances)
	for i in 1:length(instances), k in 1:length(solver_params[:algorithm])
		for t in 1:solver_params[:total_threads][k]
			params = Dict()
			params[:read_instance] = :read_bospp_xavier
			params[:read_instance_input] = instances[i]
			params[:problem_type] = "BOSPP"
			params[:instance_type] = "Xavier"
			params[:instance_id] = instances[i]
			params[:algorithm] = solver_params[:algorithm][k]
			params[:type_of_operation] = solver_params[:type_of_operation][k]
			params[:mip_emphasis] = solver_params[:mip_emphasis][k]
			params[:type_of_boost] = solver_params[:type_of_boost][k]
			params[:turn_off_cplex_heur] = solver_params[:turn_off_cplex_heur][k]
			params[:timelimit] = Inf
			params[:total_threads] = solver_params[:total_threads][k]
			params[:thread_id] = t
			push!(queue_, params)
		end
	end

	for j in 1:20, i in 1:2, k in 1:length(solver_params[:algorithm])
		for t in 1:solver_params[:total_threads][k]
			params = Dict()
			if i == 1
				params[:read_instance] = :read_bokp_hadi
				params[:problem_type] = "BOKP"
			else
				params[:read_instance] = :read_boap_hadi
				params[:problem_type] = "BOAP"
			end
			params[:read_instance_input] = j
			params[:instance_type] = "Hadi"
			params[:instance_id] = "$(j)dat.txt"
			params[:algorithm] = solver_params[:algorithm][k]
			params[:type_of_operation] = solver_params[:type_of_operation][k]
			params[:mip_emphasis] = solver_params[:mip_emphasis][k]
			params[:type_of_boost] = solver_params[:type_of_boost][k]
			params[:turn_off_cplex_heur] = solver_params[:turn_off_cplex_heur][k]
			params[:timelimit] = Inf
			params[:total_threads] = solver_params[:total_threads][k]
			params[:thread_id] = t
			push!(queue_, params)
		end
	end

	try
		files = readdir("../results/")
		files = files[[contains(files[i], "Experimental_Results") for i in 1:length(files)]]
		data = readtable("../results/$(files[1])")
		for i in 2:length(files)
    		data = vcat(data, readtable("../results/$(files[i])"))
		end
		data = data[names(data)[1:21]]
		unique!(data)
		sort!(data)
		[rm("../results/$(files[i])") for i in 1:length(files)]
		writetable("../results/Experimental_Results.csv", data)
		ind = Int64[]
		for i in 1:length(queue_)
			tmp = data[data[:Problem_Type] .== queue_[i][:problem_type],:]
    		tmp = tmp[tmp[:Instance_Type] .== queue_[i][:instance_type],:]
    		tmp = tmp[tmp[:Instance_ID] .== queue_[i][:instance_id],:]
    		tmp = tmp[tmp[:Algorithm] .== queue_[i][:algorithm],:]
    		tmp = tmp[tmp[:Type_of_Operation] .== queue_[i][:type_of_operation],:]
    		tmp = tmp[tmp[:MIP_Emphasis] .== queue_[i][:mip_emphasis],:]
    		tmp = tmp[tmp[:Type_of_Boost] .== queue_[i][:type_of_boost],:]
    		tmp = tmp[tmp[:Cplex_Heur] .!= queue_[i][:turn_off_cplex_heur],:]
 			tmp = tmp[tmp[:Parallelism] .== queue_[i][:parallelism],:]
    		tmp = tmp[tmp[:Total_Threads] .== queue_[i][:total_threads],:]
    		tmp = tmp[tmp[:Thread_ID] .== queue_[i][:thread_id],:]
    		if nrow(tmp) == 0
        		push!(ind, i)
        		continue
    		end
		end
		return unique(queue_[ind])
	catch
		return unique(queue_)
	end
end

@everywhere function solve_bopinstance_exactly(params, inds::Vector{Int64})
	for i in inds
		try
			solve_bopinstance_exactly(params[i])
		catch
		end
	end
end

@everywhere function generate_parallel_queues(inds::Vector{Int64})
	procs_::Vector{Int64} = procs()
	inds_ = Dict()
	i::Int64 = 1
	count::Int64 = 1
    while i <= length(inds)
    	p = procs_[count]
    	if p != myid() || length(procs_) == 1
        	if p in keys(inds_)
        		push!(inds_[p], inds[i])
        	else
        		inds_[p] = [inds[i]]
        	end
        	i += 1
       	end
       	count += 1
       	if count > length(procs_)
       		count = 1
       	end
    end
	inds_
end

@everywhere function run_parallel_experiments_on_all_instances()
	procs_::Vector{Int64} = procs()
	@sync begin
        for p in procs_
            @async begin
            	remotecall_fetch(testing_all_methods, p, instance, true_non_dom_sols)
            end
        end
    end
    queue = generate_queue_for_detailed_experiments_on_all_instances()
    inds = generate_parallel_queues([1:length(queue)...])
	@sync begin
        for p in procs_
        	if p != myid() || length(procs_) == 1
            	@async begin
            		remotecall_fetch(solve_bopinstance_exactly, p, queue, inds[p])
            	end
            end
        end
    end
end

#####################################################################
# Execution                                                         #
#####################################################################

#run_parallel_experiments_on_all_instances()


text: Base.:(==)(a1::VLAgentModel, a2::VLAgentModel) = (a1.agent_id == a2.agent_id)

ext: jl


text:
using Images, ImageView

function main()
    img = load("julia.jpeg")
    grayimg=Gray.(img)
    save("gray.png",grayimg)

    img = rand(368, 224, 3)
    # HxWxC ???
    #No! Actual this generates 3-images!!!
    save("img368x224.png",img)

    #img = rand(3,368,224)
    #save("img_3x368x224.png",img) dangerous do not do this

    color_g = 0.7 * ones(3,224,448) #CxHxW
    color_g = colorview(RGB, color_g)
    save("color_g.png",color_g)

    ch=color_g[120,400]
    #@show ch
    #@show red(ch)
    #@show blue(ch)
    #@show green(ch)
    #@show colorview(RGB,Gray.(rand(128,224)),Gray.(rand(128,224)),Gray.(rand(128,224)))
    img=colorview(RGB,rand(3,128,224))
    save("randgen.png",img)
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# arbres binaires de recherche mutables

mutable struct BST
    left
    value
    right
end

# principe : on modifie l'arbre *et* on renvoie toujours sa racine
#            ainsi on traite agréablement le cas de l'arbre vide
function add(x, t)
    if t == nothing
	return BST(nothing, x, nothing)
    end
    if x < t.value
	t.left = add(x, t.left)
    elseif x > t.value
	t.right = add(x, t.right)
    end
    return t
end

function mem(x, a)
    if x == a.value
        true
    elseif x < a.value && a.left != nothing
        mem(x, a.left)
    elseif a.right != nothing
        mem(x, a.right)
    else
        false
    end
end

function print_tree(a)
    if a == nothing
        return
    end
    print("(")
    if a.left != nothing
        print_tree(a.left)
    end
    print(a.value)
    if a.right != nothing
        print_tree(a.right)
    end
    print(")")
end

function main()
    tree = nothing
    for i = 1 : 9
        x = (55 * i) % 34
        tree = add(x, tree)
        print_tree(tree)
        print("\n")
    end
    if mem(8, tree) && !mem(0, tree) && mem(32, tree) && !mem(22, tree)
        println("ok")
    end
    tree = add(42, tree)
    tree = add(-1, tree)
    print_tree(tree)
    print("\n")
end

main()

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
@testset "Indexing" begin
    @testset "Linear getindex() on SVector" begin
        sv = SVector{4}(4,5,6,7)

        # Tuple
        @test (@inferred getindex(sv, (4,3,2,1))) === SVector((7,6,5,4))

        # Colon
        @test (@inferred getindex(sv,:)) === sv
    end

    @testset "Linear getindex()/setindex!() on MVector" begin
        vec = @SVector [4,5,6,7]

        # Tuple
        mv = MVector{4,Int}()
        @test (mv[(1,2,3,4)] = vec; (@inferred getindex(mv, (4,3,2,1)))::MVector{4,Int} == MVector((7,6,5,4)))

        # Colon
        mv = MVector{4,Int}()
        @test (mv[:] = vec; (@inferred getindex(mv, :))::MVector{4,Int} == MVector((4,5,6,7)))
    end

    @testset "Linear getindex()/setindex!() with a tuple on an Array" begin
        v = [11,12,13]
        m = [1.0 2.0; 3.0 4.0]

        @test v[(2,3)] === SVector(12, 13)
        @test m[(2,3)] === SVector(3.0, 2.0)

        @test (v[(2,3)] = [22,23]; (v[2] == 22) & (v[3] == 23))

    end

    @testset "2D getindex() on SMatrix" begin
        sm = @SMatrix [1 3; 2 4]

        # Scalar
        @test sm[1,1] === 1
        @test sm[2,1] === 2
        @test sm[1,2] === 3
        @test sm[2,2] === 4

        # Tuple, scalar
        @test (@inferred getindex(sm, (2,1), (2,1))) === @SMatrix [4 2; 3 1]
        @test (@inferred getindex(sm, 1, (1,2))) === @SVector [1,3]
        @test (@inferred getindex(sm, (1,2), 1)) === @SVector [1,2]

        # Colon
        @test (@inferred getindex(sm, :, :)) === @SMatrix [1 3; 2 4]
        @test (@inferred getindex(sm, (2,1), :)) === @SMatrix [2 4; 1 3]
        @test (@inferred getindex(sm, :, (2,1))) === @SMatrix [3 1; 4 2]
        @test (@inferred getindex(sm, 1, :)) === @SVector [1,3]
        @test (@inferred getindex(sm, :, 1)) === @SVector [1,2]
    end

    @testset "2D getindex()/setindex! on MMatrix" begin
        sm = @MMatrix [1 3; 2 4]

        # Tuple, scalar
        @test (mm = MMatrix{2,2,Int}(); mm[(2,1),(2,1)] = sm[(2,1),(2,1)]; (@inferred getindex(mm, (2,1), (2,1)))::MMatrix == @MMatrix [4 2; 3 1])
        @test (mm = MMatrix{2,2,Int}(); mm[1,(1,2)] = sm[1,(1,2)]; (@inferred getindex(mm, 1, (1,2)))::MVector == @MVector [1,3])
        @test (mm = MMatrix{2,2,Int}(); mm[(1,2),1] = sm[(1,2),1]; (@inferred getindex(mm, (1,2), 1))::MVector == @MVector [1,2])

        # Colon
        @test (mm = MMatrix{2,2,Int}(); mm[:,:] = sm[:,:]; (@inferred getindex(mm, :, :))::MMatrix == @MMatrix [1 3; 2 4])
        @test (mm = MMatrix{2,2,Int}(); mm[(2,1),:] = sm[(2,1),:]; (@inferred getindex(mm, (2,1), :))::MMatrix == @MMatrix [2 4; 1 3])
        @test (mm = MMatrix{2,2,Int}(); mm[:,(2,1)] = sm[:,(2,1)]; (@inferred getindex(mm, :, (2,1)))::MMatrix == @MMatrix [3 1; 4 2])
        @test (mm = MMatrix{2,2,Int}(); mm[1,:] = sm[1,:]; (@inferred getindex(mm, 1, :))::MVector == @MVector [1,3])
        @test (mm = MMatrix{2,2,Int}(); mm[:,1] = sm[:,1]; (@inferred getindex(mm, :, 1))::MVector == @MVector [1,2])
    end

    @testset "2D getindex() with tuples on an Array" begin
        m = [1.0 2.0; 3.0 4.0]

        @test m[(1,2), (1,2)] === @SMatrix [1.0 2.0; 3.0 4.0]
        @test m[1, (1,2)] === @SVector [1.0, 2.0]
        @test m[(1,2), 1] === @SVector [1.0, 3.0]
    end

    @testset "3D scalar indexing" begin
        sa = SArray{(2,2,2), Int}([i*j*k for i = 1:2, j = 2:3, k=3:4])

        @test sa[1,1,2] === 8
        @test sa[1,2,1] === 9
        @test sa[2,1,1] === 12

        ma = MArray{(2,2,2), Int}()
        @test (ma[1,1,2] = 8; ma[1,1,2] === 8)
        @test (ma[1,2,1] = 9; ma[1,2,1] === 9)
        @test (ma[2,1,1] = 12; ma[2,1,1] === 12)
    end

    @testset "4D scalar indexing" begin
        sa = SArray{(2,2,2,2), Int}([i*j*k*l for i = 1:2, j = 2:3, k=3:4, l=4:5])

        @test sa[1,1,1,2] === 30
        @test sa[1,1,2,1] === 32
        @test sa[1,2,1,1] === 36
        @test sa[2,1,1,1] === 48

        ma = MArray{(2,2,2,2), Int}()
        @test (ma[1,1,1,2] = 30; ma[1,1,1,2] === 30)
        @test (ma[1,1,2,1] = 32; ma[1,1,2,1] === 32)
        @test (ma[1,2,1,1] = 36; ma[1,2,1,1] === 36)
        @test (ma[2,1,1,1] = 48; ma[2,1,1,1] === 48)
    end

end


text:
"""
    create_multi_year_network_data(case, number_of_hours, number_of_scenarios, number_of_years)

    Using the input case (case6, cigre, ...) create multi-year network data using the add_dimensions!(...).

    Dimension hierarchy is:
    year{...
        scenario{...
            hour{...}
        }
    }

"""
function create_multi_year_network_data(case, number_of_hours, number_of_scenarios, number_of_years; kwargs...)
    my_data = Dict{String, Any}("multinetwork"=>true, "name"=>case, "nw"=>Dict{String,Any}(), "per_unit"=>true)
    if case == "case6"
        base_file = normpath(@__DIR__,"..","..","test","data","case6","case6_")
        planning_years = [2030, 2040, 2050]

        _FP.add_dimension!(my_data, :hour, number_of_hours)
        scenario = Dict(s => Dict{String,Any}("probability"=>1/number_of_scenarios) for s in 1:number_of_scenarios)
        _FP.add_dimension!(my_data, :scenario, scenario, metadata = Dict{String,Any}("mc"=>get(kwargs, :mc, true)))
        _FP.add_dimension!(my_data, :year, number_of_years; metadata = Dict{String,Any}("scale_factor"=>10))
    elseif case == "case67"
        base_file = normpath(@__DIR__,"..","..","test","data","case67","case67_tnep_")
        planning_years = [2030, 2040, 2050]

        _FP.add_dimension!(my_data, :hour, number_of_hours)
        data_years = [2017, 2018, 2019]
        start = [1483228800000, 1514764800000, 1546300800000]
        scenario = Dict(s => Dict{String,Any}(
                "probability"=>1/number_of_scenarios,
                "year" => data_years[s],
                "start" => start[s]
            ) for s in 1:number_of_scenarios)
        _FP.add_dimension!(my_data, :scenario, scenario)
        _FP.add_dimension!(my_data, :year, number_of_years; metadata = Dict{String,Any}("scale_factor"=>10))
    else
        error("Case \"$(case)\" not (yet) supported.")
    end

    cost_scale_factor = get(kwargs, :cost_scale_factor, 1.0)
    for year_idx = 1 : number_of_years
        year = planning_years[year_idx]
        file = base_file * "$year" * ".m"
        data = _FP.parse_file(file)
        data["dim"] = _FP.dim(my_data)
        _FP.scale_data!(data; year_idx, cost_scale_factor)
        add_one_year!(my_data, case, data, year_idx)
    end

    return my_data
end

function add_one_year!(my_data, case, data, year_idx)
    number_of_nws = _FP.dim_length(data, :hour) * _FP.dim_length(data, :scenario)
    nw_id_offset = number_of_nws * (year_idx - 1)
    if case == "case6"
        data, loadprofile, genprofile = create_profile_data_italy!(data)
    elseif case == "case67"
        data, loadprofile, genprofile = create_profile_data_germany!(data)
    else
        error("Case \"$(case)\" not (yet) supported.")
    end
    time_series = create_profile_data(number_of_nws, data, loadprofile, genprofile)
    mn_data = _FP.make_multinetwork(data, time_series; number_of_nws, nw_id_offset)
    _FP.import_nws!(my_data, mn_data)
    return my_data
end


text:
mutable struct DecisionTree <: Classifier
    root      :: Union{Nothing, Node}
    topredict
    maxdepth  :: UInt
    criterion
    attrfilter
end

DecisionTree(;maxdepth = typemax(Int), criterion = gini, attrfilter = identity) = begin
    DecisionTree(nothing, nothing, maxdepth, criterion, attrfilter)
end

fit!(dt::DecisionTree, df::DataFrame, topredict) = begin
    dt.topredict = topredict
    dt.root = Node(df, setdiff(propertynames(df), [topredict]), 0, dt)
    dt
end

predict(dt::DecisionTree, df) = predict(dt.root, df)


text:
using FastqSubsample

ifastq = ARGS[1]
ofastq = ARGS[2]
nSubsample = parse(Int, ARGS[3])
seed = parse(Int, ARGS[4])
fastqsubsample(ifastq, ofastq, nSubsample, seed = seed)

max_stars_repo_name: yuifu/FastqSubsample

max_issues_repo_name: yuifu/FastqSubsample

max_forks_repo_name: yuifu/FastqSubsample


text:
function test_constructW(D, T)
    W = VoiceConversion.constructW(D, T)

    @test issparse(W)
    @test size(W) == (2D*T, D*T)

    for t=1:T
        # coef. for static feature
        s = 2D*(t-1)+1
        e = s + D - 1
        @test W[s:e,(t-1)*D+1:t*D] == speye(D)
        for i=1:T
            if i != t
                @test W[s:e, (i-1)*D+1:i*D] == spzeros(D, D)
            end
        end

        # coef. for dynamic feature
        s = s + D
        e = e + D
        if t >= 2
            @test W[s:e, (t-2)*D+1:(t-1)*D] == -0.5*speye(D)
        end
        if t < T
            @test W[s:e, t*D+1:(t+1)*D] == 0.5*speye(D)
        end

        for i=1:T
            if i != t-1 && i != t+1
                @test W[s:e, (i-1)*D+1:i*D] == spzeros(D, D)
            end
        end
    end
end

function test_trajectory(T)
    modelpath = joinpath(Pkg.dir("VoiceConversion"), "test", "models",
                         "clb_to_slt_gmm32_order40_diff_with_delta.jld")
    gmm = load(modelpath)
    @assert gmm["diff"]

    mapper = GMMMap(gmm["weights"], gmm["means"], gmm["covars"])
    mapper = TrajectoryGMMMap(mapper, T)
    @test length(mapper) == T

    D = div(size(gmm["means"], 1), 4)
    @test dim(mapper) == 2D # must contains delta
    @test ncomponents(mapper) == length(gmm["weights"])
    @test size(mapper) == (2D, T)
end

println("testing: GMM-based trajectory parameter conversion")
test_constructW(30, 40)
test_trajectory(100)


text:
@testset "Testing lqcontrol.jl" begin

    rough_kwargs = Dict(:atol => 1e-13, :rtol => 1e-4)

    # set up
    q    = 1.
    r    = 1.
    rf   = 1.
    a    = .95
    b    = -1.
    c    = .05
    β    = .95
    n    = 0.
    capT = 1
    lq_scalar = LQ(q, r, a, b, c, n, β, capT, rf)

    Q  = [0. 0.; 0. 1]
    R  = [1. 0.; 0. 0]
    rf = eye(2) .* 100
    A  = fill(0.95, 2, 2)
    B  = fill(-1.0, 2, 2)
    lq_mat = LQ(Q, R, A, B, bet=β, capT=capT, rf=rf)

    @testset "Test scalar sequences with exact by hand solution" begin
        x0 = 2.0
        x_seq, u_seq, w_seq = compute_sequence(lq_scalar, x0)
        # solve by hand
        u_0 = (-2.*lq_scalar.A*lq_scalar.B*lq_scalar.bet*lq_scalar.rf) /
           (2.*lq_scalar.Q+lq_scalar.bet*lq_scalar.rf*2lq_scalar.B^2)*x0
        x_1 = lq_scalar.A * x0 + lq_scalar.B * u_0 + w_seq[end]

        @test isapprox(u_0[1], u_seq[end]; rough_kwargs...)
        @test isapprox(x_1[1], x_seq[end]; rough_kwargs...)
    end

    @testset "test matrix solutions" begin
        x0 = randn(2) * 25
        x_seq, u_seq, w_seq = compute_sequence(lq_mat, x0)

        @test isapprox(sum(u_seq), 0.95 * sum(x0); rough_kwargs...)
        @test isapprox(x_seq[:,end], zeros(x0); rough_kwargs...)
    end

    @testset "test stationary matrix" begin
        x0 = randn(2) .* 25
        P, F, d = stationary_values(lq_mat)

        f_answer = [-.95 -.95; 0. 0.]
        p_answer = [1. 0; 0. 0.]

        val_func_lq = (x0' * P * x0)[1]
        val_func_answer = x0[1]^2

        @test isapprox(f_answer, F; rough_kwargs...)
        @test isapprox(val_func_lq, val_func_answer; rough_kwargs...)
        @test isapprox(p_answer, P; rough_kwargs...)
    end

    @testset "test runs a (n,k,j) = (2,1,1) model" begin
        # == Model parameters == #
        r = 0.05
        bet = 1 / (1 + r)
        t = 45
        c_bar = 2.0
        sigma = 0.25
        mu = 1.0
        q = 1e6

        # == Formulate as an LQ problem == #
        Q = 1.0
        R = zeros(2, 2)
        Rf = zeros(2, 2); Rf[1, 1] = q
        A = [1.0+r -c_bar+mu;
             0.0 1.0]
        B = [-1.0, 0.0]
        C = [sigma, 0.0]

        # == Compute solutions and simulate == #
        lq = LQ(Q, R, A, B, C; bet=bet, capT=t, rf=Rf)
        x0 = [0.0, 1.0]
        xp, up, wp = compute_sequence(lq, x0)
        @test true == true  # just assert true if we made it to this point
    end

end  # @testset


text:
using Documenter, PowerSystemCaseBuilder
# import DataStructures: OrderedDict
# using Literate

if isfile("docs/src/howto/.DS_Store.md")
    rm("docs/src/howto/.DS_Store.md")
end

makedocs(
    sitename = "PowerSystemCaseBuilder.jl",
    format = Documenter.HTML(
        mathengine = Documenter.MathJax(),
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    modules = [PowerSystemCaseBuilder],
    strict = true,
    authors = "Sourabh Dalvi",
    pages = Any[
        "Introduction" => "index.md",
    ],
)

deploydocs(
    repo = "github.com/NREL-SIIP/PowerSystemCaseBuilder.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#"],
)

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
__precompile__()

module GadgetIO

    using Printf

    function output_time(t1, t2)
        return @sprintf("%0.3e", Float64((t2-t1))*1.e-9)
    end

    """
        AbstractGadgetHeader

    Abstract type to unify snapshot and subfind headers.
    """
    abstract type AbstractGadgetHeader end

    # functions to read snapshots
    include(joinpath("read_snapshot", "utility", "domain_slices.jl"))
    include(joinpath("read_snapshot", "read_header.jl"))
    include(joinpath("read_snapshot", "read_info.jl"))
    include(joinpath("read_snapshot", "utility", "allocate_dict.jl"))
    include(joinpath("read_snapshot", "utility", "snapshot_utilities.jl"))
    include(joinpath("read_snapshot", "read_format_1.jl"))
    include(joinpath("read_snapshot", "read_format_2.jl"))
    include(joinpath("read_snapshot", "read_snapshot.jl"))
    include(joinpath("read_subfind",  "read_subfind.jl"))
    include(joinpath("read_subfind",  "read_halo_prop.jl"))
    include(joinpath("read_subfind",  "filter_subfind.jl"))
    include(joinpath("read_snapshot", "particles_in_box",  "read_key_files.jl"))
    include(joinpath("read_snapshot", "particles_in_box",  "peano_hilbert.jl"))
    include(joinpath("read_snapshot", "particles_in_box",  "utility.jl"))
    include(joinpath("read_snapshot", "particles_in_box",  "read_particle_in_box.jl"))
    include(joinpath("read_snapshot", "particles_in_geometry", "abstract_geometry.jl"))
    include(joinpath("read_snapshot", "particles_in_geometry", "cylinder.jl"))
    include(joinpath("read_snapshot", "particles_in_geometry", "sphere.jl"))
    include(joinpath("read_snapshot", "particles_in_geometry", "cube.jl"))
    include(joinpath("read_snapshot", "particles_in_geometry", "read_particles_in_geometry.jl"))
    include(joinpath("read_snapshot", "particles_in_halo", "read_particles_by_id.jl"))
    include(joinpath("read_snapshot", "particles_in_halo", "read_particles_in_halo.jl"))
    include(joinpath("read_snapshot", "distributed_files", "find_read_positions.jl"))
    include(joinpath("read_snapshot", "distributed_files", "read_distributed_files.jl"))

    # functions to write snapshots
    include(joinpath(dirname(@__FILE__), "write_snapshot", "write_snap.jl"))

    export  AbstractGadgetHeader,
            SnapshotHeader, SubfindHeader,
            InfoLine,       # types
            head_to_dict,
            snap_to_dict,
            head_to_obj,
            print_blocks,
            read_info,
            block_present,
            read_snap,
            read_block,      # similar to readnew.pro by Klaus Dolag
            read_header,

            # large simulations
            read_particles_in_box,
            read_particles_in_volume,
            read_particles_in_geometry,
            read_particles_in_halo,
            get_index_list,
            get_npart_to_read,
            filter_ids,
            read_blocks_over_all_files,
            find_read_positions,
            save_read_positions,
            load_read_positions,

            # geometries
            AbstractGadgetGeometry,
            GadgetCylinder,
            GadgetSphere,
            GadgetCube,

            # domain slices
            filter_cube,
            filter_cylinder,

            # subfind read
            HaloID,
            SubfindHeader,
            read_subfind_header,
            read_subfind,
            find_most_massive_halo,
            filter_subfind,
            read_ids_in_halo,
            read_particles_by_id,
            read_halo_prop,
            read_halo_prop_and_id,

            # write snapshot functions
            write_header,
            write_block,
            write_info_block

end # module


text:
export ad, mad, mad_core

"""
`mad_model(G)` returns the solved linear program whose optimum
value is `mad(G)`.
"""
function mad_model(G::SimpleGraph)
  EE = elist(G)
  VV = vlist(G)

  n = length(VV)
  m = length(EE)

  MOD = Model(get_solver())

  # variables
  @variable(MOD, z >= 0)
  @variable(MOD, x[i=1:n,j=1:m; in(VV[i],EE[j])] >= 0 )

  # vertex constraints
  for i=1:n
    v = VV[i]
    # get all edges that have v as an end point
    j_iterator = (j for j=1:m if EE[j][1]==v || EE[j][2]==v)
    j_list = collect(j_iterator)
    @constraint(MOD, sum(x[i,j] for j in j_list) <= z)
  end

  # edge constraints
  for j=1:m
    v,w = EE[j]
    i1 = first(findall( x -> x==v, VV))
    i2 = first(findall( x -> x==w, VV))
    @constraint(MOD, x[i1,j]+x[i2,j] == 2)
  end

  @objective(MOD, Min, z)
  optimize!(MOD)
  return MOD, x, VV, EE
end

"""
`mad(G)` computes the maximum average degree of `G`.
See also `mad_core`.
"""
function mad(G::SimpleGraph)
  if cache_check(G,:mad)
    return cache_recall(G,:mad)
  end
  MOD, x, VV, EE = mad_model(G)
  mval = objective_value(MOD)
  cache_save(G,:mad,mval)
  return mval
end

"""
`ad(G)` is the average degree of a vertex in `G`.
"""
ad(G::SimpleGraph) = 2*NE(G)/NV(G)

"""
`mad_core(G)` returns a subgraph `H` of `G` for which
`ad(H)==mad(G)`.
"""
function mad_core(G::SimpleGraph)
  if cache_check(G,:mad_core)
    return cache_recall(G,:mad_core)
  end
  T = vertex_type(G)
  GG = deepcopy(G)
  while true
    if NV(GG)==0
      error("Bad graph or something went wrong")
    end
    # solve the LP for GG
    MOD, x, VV, EE = mad_model(GG)
    n = length(VV)
    m = length(EE)
    err = 0.1/n
    mval = objective_value(MOD)

    # if balanced, return
    if abs(ad(GG)-mval) <= err
      cache_save(G,:mad_core,GG)
      return GG
    end

    # otherwise, find and destroy defective vertices
    for i=1:n
      v = VV[i]
      total = 0.0
      for j=1:m
        e = EE[j]
        if in(v,e)
          val = value(x[i,j])
          total += val
        end
      end
      if total < mval - err   # this is a defective vertex
        delete!(GG,v)
      end
    end
  end # end while
  error("This can't happen (but it did)")
  cache_save(G,:mad_core,GG)
  return GG
end


text:
using ImageBase
using Statistics
using Test

@testset "Reductions" begin
    _abs(x::Colorant) = mapreducec(abs, +, 0, x)

    A = rand(5,5,3)
    img = colorview(RGB, PermutedDimsArray(A, (3,1,2)))
    s12 = sum(img, dims=(1,2))
    @test eltype(s12) <: RGB
    A = [NaN, 1, 2, 3]
    @test meanfinite(A, 1) ≈ [2]
    A = [NaN 1 2 3;
         NaN 6 5 4]
    mf = meanfinite(A, 1)
    @test isnan(mf[1])
    @test mf[1,2:end] ≈ [3.5,3.5,3.5]
    @test meanfinite(A, 2) ≈ reshape([2, 5], 2, 1)
    @test meanfinite(A, (1,2)) ≈ [3.5]
    @test minfinite(A) == 1
    @test maxfinite(A) == 6
    @test maxabsfinite(A) == 6
    A = rand(10:20, 5, 5)
    @test minfinite(A) == minimum(A)
    @test maxfinite(A) == maximum(A)
    A = reinterpret(N0f8, rand(0x00:0xff, 5, 5))
    @test minfinite(A) == minimum(A)
    @test maxfinite(A) == maximum(A)
    A = rand(Float32,3,5,5)
    img = colorview(RGB, A)
    dc = meanfinite(img, 1)-reshape(reinterpretc(RGB{Float32}, mean(A, dims=2)), (1,5))
    @test maximum(map(_abs, dc)) < 1e-6
    dc = minfinite(img)-RGB{Float32}(minimum(A, dims=(2,3))...)
    @test _abs(dc) < 1e-6
    dc = maxfinite(img)-RGB{Float32}(maximum(A, dims=(2,3))...)
    @test _abs(dc) < 1e-6
end


text:
# ## Getting started

using MLJ, RDatasets, PrettyPrinting, Random

# We start by generating a 2D cloud of points

Random.seed!(3203)
X = randn(20, 2)
y = vcat(-ones(10), ones(10))

X = MLJ.table(X)
y = categorical(y);

# Let's fit a SVM classifier

@load SVC pkg=LIBSVM

svc_mdl = SVC()
svc = machine(svc_mdl, X, y)

fit!(svc);

# As usual we can check how it performs

ypred = predict(svc, X)
misclassification_rate(ypred, y)

# Not bad.

# ### Basic tuning

# As usual we could tune the model, for instance the penalty encoding the tradeoff between margin width and misclassification:

rc = range(svc_mdl, :cost, lower=0.1, upper=5)
tm = TunedModel(model=svc_mdl, ranges=[rc], tuning=Grid(resolution=10),
                resampling=CV(nfolds=3, rng=33), measure=misclassification_rate)
mtm = machine(tm, X, y)

fit!(mtm)

ypred = predict(mtm, X)
misclassification_rate(ypred, y)

# You could also change the kernel etc.

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
@testset "Tranformations" begin

@testset "Rotations" begin

end

end


text:
export latex

using REPL
using Printf

const _SUPERSCRIPTNUMS = collect("⁰¹²³⁴⁵⁶⁷⁸⁹")
subscript(i::Integer) = (i<0 ? "₋" : "") * join('₀'+d for d in reverse(digits(abs(i))))
superscript(i::Integer) = (i<0 ? "⁻" : "") * join(_SUPERSCRIPTNUMS[d+1] for d in reverse(digits(abs(i))))

function printspaced(io::IO,v,trailing_space=false)
    isempty(v) && return
    print(io,v[1])
    iprinted = 1
    exponent = 1
    for i = 2:length(v)
        if v[i]==v[iprinted]
            exponent += 1
            # if the last element is repeated, print the exponent
            i == length(v) && print(io,superscript(exponent))
        else
            # print exponent for previous object
            exponent>1 && print(io,superscript(exponent))
            exponent = 1
            print(io," ")
            print(io,v[i])
            iprinted = i
        end
    end
    trailing_space && print(io," ")
end

Base.print(io::IO,A::BaseOpProduct) = printspaced(io,A.v)
function Base.print(io::IO,ii::QuIndex)
    if isnoindex(ii)
        return
    elseif isintindex(ii)
        print(io,ii.num)
    else
        print(io,ii.sym)
        ii.num != typemin(ii.num) && print(io,subscript(ii.num))
    end
end

_print(io::IO,A::BaseOperator) = print(io,A.name,BaseOpType_sym[Int(A.t)],"(",A.inds...,")")
Base.print(io::IO,A::BaseOperator) = _print(io, unalias(A))
Base.print(io::IO,A::δ) = print(io,"δ(",A.iA,A.iB,")")
Base.print(io::IO,A::Param) = print(io, A.name, A.state=='c' ? "*" : "", length(A.inds)==0 ? "" : "($(A.inds...))")
Base.print(io::IO,A::ExpVal) = print(io,"⟨", A.ops, "⟩")
Base.print(io::IO,A::Corr) = print(io,"⟨", A.ops, "⟩c")

function Base.print(io::IO,A::QuTerm)
    if A.nsuminds > 0
        print(io,"∑")
        for ii=1:A.nsuminds
            print(io,subscript(ii))
        end
        print(io," ")
    end
    isempty(A.δs) || printspaced(io,A.δs,true)
    isempty(A.params) || printspaced(io,A.params,true)
    isempty(A.expvals) || printspaced(io,A.expvals,true)
    isempty(A.corrs) || printspaced(io,A.corrs,true)
    isempty(A.bares) || print(io,A.bares)
end

function print_as_connector(io,str::String,print_connector::Bool)
    if print_connector
        if !isempty(str) && str[1] == '-'
            print(io," - ",str[2:end])
        else
            print(io," + ",str)
        end
    else
        print(io,str)
    end
end

function print_term_scalar(io::IO,t::QuTerm,s::Number, print_connector::Bool)
    if isempty(t)
        print_as_connector(io,string(s),print_connector)
        return
    elseif s == -one(s)
        str = "-"
    elseif s == 1
        str = ""
    else
        str = numstring(s)*" "
    end
    print_as_connector(io,str,print_connector)
    print(io,t)
end

numstring(x::Number) = @sprintf "%g" x
numstring(x::Rational) = denominator(x)==1 ? string(numerator(x)) : string(x)
numstring(x::Complex) = ((r,i) = reim(x); iszero(i) ? numstring(r) : (iszero(r) ? numstring(i)*"i" : "($(numstring(r))$(i<0 ? '-' : '+')$(numstring(abs(i)))i)"))

function Base.print(io::IO,A::QuExpr)
    if isempty(A.terms)
        print(io,0)
    else
        for (ii,(t,s)) in enumerate(sort!(Tuple.(collect(A.terms))))
            print_connector = ii != 1
            print_term_scalar(io,t,s, print_connector)
        end
    end
end

Base.show(io::IO, A::QuantumObject) = print(io,A)
Base.show(io::IO, ::MIME"text/latex", A::QuantumObject) = (print(io,"\$"); printlatex(io,A); print(io,"\$"))

function latex(A::Union{QuantumObject,Number})
    io = IOBuffer()
    printlatex(io, A)
    String(take!(io))
end

function printlatex(io::IO,ii::QuIndex)
    if isnoindex(ii)
        return
    elseif isintindex(ii)
        print(io,ii.num)
    elseif issumindex(ii)
        print(io,"\\#_{",ii.num,"}")
    else
        print(io,ii.sym)
        ii.num != typemin(ii.num) && print(io,"_{",ii.num,"}")
    end
end

function printlatex(io::IO,v::Vector)
    isempty(v) && return
    prevA = v[1]
    exponent = 1
    for i = 2:length(v)
        if v[i]==prevA
            exponent += 1
        else
            # print previous object
            exponent > 1 && print(io,"{")
            printlatex(io,prevA)
            exponent > 1 && print(io,"}^{",exponent,"}")
            exponent = 1
            prevA = v[i]
            print(io," ")
        end
    end
    # print last object
    exponent > 1 && print(io,"{")
    printlatex(io,prevA)
    exponent > 1 && print(io,"}^{",exponent,"}")
end

latexjoin(args...) = join(latex.(args...))
latexindstr(inds::QuIndex...) = latexindstr(inds)
latexindstr(inds) = isempty(inds) ? "" : "_{$(latexjoin(inds))}"

numlatex(x::Number) = @sprintf "%g" x
numlatex(x::Rational) = denominator(x)==1 ? "$(numerator(x))" : "\\frac{$(numerator(x))}{$(denominator(x))}"
numlatex(x::Complex) = @sprintf "(%g%+gi)" real(x) imag(x)
numlatex(x::Complex{Rational{T}}) where T = @sprintf "\\left(%s%s%si\\right)" numlatex(real(x)) (imag(x)>=0 ? "+" : "-") numlatex(abs(imag(x)))

const _unicode_to_latex = Dict(v[1]=>"{$k}" for (k,v) in REPL.REPLCompletions.latex_symbols if length(v)==1)
unicode_to_latex(s::AbstractString) = join(map(c -> get(_unicode_to_latex,c,string(c)), collect(s)))
unicode_to_latex(s::QuOpName) = unicode_to_latex(string(s))

printlatex(io::IO,x::Number) = print(io, imag(x)==0 ? numlatex(real(x)) : (real(x)==0 ? numlatex(imag(x))*"i" : numlatex(x)))
_printlatex(io::IO,A::BaseOperator) = print(io,"{", unicode_to_latex(A.name), "}", latexindstr(A.inds), BaseOpType_latex[Int(A.t)])
printlatex(io::IO,A::BaseOperator) = _printlatex(io, unalias(A))
printlatex(io::IO,A::δ) = print(io, "\\delta", latexindstr(A.iA,A.iB))
printlatex(io::IO,A::Param) = print(io, unicode_to_latex(A.name), latexindstr(A.inds), A.state=='c' ? "^{*}" : "")
printlatex(io::IO,A::BaseOpProduct) = printlatex(io,A.v)
printlatex(io::IO,A::ExpVal) = (print(io, "\\langle "); printlatex(io, A.ops); print(io," \\rangle"))
printlatex(io::IO,A::Corr)   = (print(io, "\\langle "); printlatex(io, A.ops); print(io," \\rangle_{c}"))

function printlatex(io::IO,A::QuTerm)
    if A.nsuminds > 0
        print(io,"\\sum",latexindstr(sumindex.(1:A.nsuminds)))
    end
    printlatex(io,A.δs)
    printlatex(io,A.params)
    printlatex(io,A.expvals)
    printlatex(io,A.corrs)
    printlatex(io,A.bares)
end

function printlatex(io::IO,A::QuExpr)
    if isempty(A)
        print(io,0)
    else
        # convert to tuples so sorting ignores the numbers if terms are different (which they are)
        for (ii,(t,s)) in enumerate(sort!(Tuple.(collect(A.terms))))
            if isempty(t)
                printlatex(io,s)
            else
                ls = s==one(s) ? "" : (s==-one(s) ? "-" : latex(s))
                if ii > 1 && (ls=="" || (ls[1] ∉ ('-','+')))
                    print(io," + ")
                end
                print(io,ls)
                printlatex(io,t)
            end
        end
    end
end


text:
function play_trajectory(vis, problem, robot, x)
    ee_positions = get_ee_path(problem, robot, x)
    show_ee_path(vis, ee_positions)

    nₓ = robot.n_q + robot.n_v + robot.n_τ  # dimension of each mesh point
    ind_q = hcat([range(1 + (i * nₓ), length=robot.n_q) for i = (1:problem.num_knots) .- 1]...)
    q_mat = x[ind_q]

    ts = range(0, length=problem.num_knots, step=problem.dt)
    qs = [q_mat[:, i] for i = 1:size(q_mat, 2)]

    animation = MeshCat.Animation(robot.mvis, ts, qs)
    setanimation!(robot.mvis, animation)
end

function plot_results(problem::Problem, robot::Robot, x)
    nₓ = robot.n_q + robot.n_v + robot.n_τ  # dimension of each mesh point

    # Helper ranges
    ind_q = hcat([range(1 + (i * nₓ)                        , length=robot.n_q) for i = (1:problem.num_knots    ) .- 1]...)
    ind_v = hcat([range(1 + (i * nₓ) + robot.n_q            , length=robot.n_v) for i = (1:problem.num_knots    ) .- 1]...)
    ind_τ = hcat([range(1 + (i * nₓ) + robot.n_q + robot.n_v, length=robot.n_τ) for i = (1:problem.num_knots - 1) .- 1]...)

    ts = range(0, length=problem.num_knots, step=problem.dt)

    plts = plot(ts         , x[ind_q]', lw=2, legend=nothing, xlabel="time [s]", ylabel="position [rad]"),
           plot(ts         , x[ind_v]', lw=2, legend=nothing, xlabel="time [s]", ylabel="velocity [rad/s]"),
           plot(ts[1:end-1], x[ind_τ]', lw=2, legend=nothing, xlabel="time [s]", ylabel="torque [Nm]")

    plot(
        plts...,
        layout=grid(3,1),
        size=(600,3*400),
        left_margin=20px,
        bottom_margin=20px,
    )
end

function plot_log(solver_log::SolverLog)

    plt_1 = plot(
        title="feasibility", title_location=:left,
        xlabel="iterations",
        ylabel="error",
        yaxis=(:log10, [1e-11, 1e4]),
        legend=nothing,
    )

    plot!(solver_log.abs_feas_error, c=1, lw=2)
    scatter!([length(solver_log.abs_feas_error)], solver_log.abs_feas_error[end:end], c=1, marker=:star4)

    plt_2 = plot(
        title="objective", title_location=:left,
        xlabel="iterations",
        ylabel="cost (1e4)",
        legend=nothing,
    )

    obj_value = solver_log.obj_value / 1e4

    plot!(obj_value, c=1, lw=2)
    scatter!([length(obj_value)], obj_value[end:end], c=1, marker=:star4)

    plot(
        plt_1, plt_2,
        size=(600,300),
        # left_margin=20px,

        fontfamily = "Times",
        guidefontsize = 12,
        legendfontsize = 12,
        tickfontsize = 10,
        titlefontsize = 15,
    )
end

export
    play_trajectory,
    plot_results,
    plot_log


text:
struct JoinInner    # inner
end
struct JoinBoth     # inner
end
struct JoinAll      # outer
end
struct JoinOuter    # outer
end
struct JoinLeft     # left
end
struct JoinRight    # right
end

"""
# Joins/Column-binding

`TS` objects can be combined together column-wise using `Index` as the
column key. There are four kinds of column-binding operations possible
as of now. Each join operation works by performing a Set operation on
the `Index` column and then merging the datasets based on the output
from the Set operation. Each operation changes column names in the
final object automatically if the operation encounters duplicate
column names amongst the TS objects.

The following join types are supported:

`join(ts1::TS, ts2::TS, ::Type{JoinInner})` and
`join(ts1::TS, ts2::TS, ::Type{JoinBoth})`

a.k.a. inner join, takes the intersection of the indexes of `ts1` and
`ts2`, and then merges the columns of both the objects. The resulting
object will only contain rows which are present in both the objects'
indexes. The function will rename columns in the final object if
they had same names in the TS objects.

`join(ts1::TS, ts2::TS, ::Type{JoinOuter})` and
`join(ts1::TS, ts2::TS, ::Type{JoinAll})`:

a.k.a. outer join, takes the union of the indexes of `ts1` and `ts2`
before merging the other columns of input objects. The output will
contain rows which are present in all the input objects while
inserting `missing` values where a row was not present in any of the
objects. This is the default behaviour if no `JoinType` object is
provided.

`join(ts1::TS, ts2::TS, ::Type{JoinLeft})`:

Left join takes the index values which are present in the left
object `ts1` and finds matching index values in the right object
`ts2`. The resulting object includes all the rows from the left
object, the column values from the left object, and the values
associated with matching index rows on the right. The operation
inserts `missing` values where in the unmatched rows of the right
object.

`join(ts1::TS, ts2::TS, ::Type{JoinRight})`

Right join, similar to left join but works in the opposite
direction. The final object contains all the rows from the right
object while inserting `missing` values in rows missing from the left
object.

The default behaviour is to assume `JoinAll` if no `JoinType` object
is provided to the `join` method.

`cbind` is an alias for `join` method.

# Examples
```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random;

julia> random(x) = rand(MersenneTwister(123), x);

julia> dates = collect(Date(2017,1,1):Day(1):Date(2017,1,10));

julia> ts1 = TS(random(length(dates)), dates)
julia> show(ts1)
(10 x 1) TS with Dates.Date Index

 Index       x1
 Date        Float64
───────────────────────
 2017-01-01  0.768448
 2017-01-02  0.940515
 2017-01-03  0.673959
 2017-01-04  0.395453
 2017-01-05  0.313244
 2017-01-06  0.662555
 2017-01-07  0.586022
 2017-01-08  0.0521332
 2017-01-09  0.26864
 2017-01-10  0.108871

julia> dates = collect(Date(2017,1,1):Day(1):Date(2017,1,30));

julia> ts2 = TS(random(length(dates)), dates);
julia> show(ts2)
(30 x 1) TS with Dates.Date Index

 Index       x1
 Date        Float64
───────────────────────
 2017-01-01  0.768448
 2017-01-02  0.940515
 2017-01-03  0.673959
 2017-01-04  0.395453
 2017-01-05  0.313244
 2017-01-06  0.662555
 2017-01-07  0.586022
 2017-01-08  0.0521332
     ⋮           ⋮
 2017-01-23  0.281066
 2017-01-24  0.792931
 2017-01-25  0.20923
 2017-01-26  0.918165
 2017-01-27  0.614255
 2017-01-28  0.802665
 2017-01-29  0.555668
 2017-01-30  0.940782
        14 rows omitted

# join on all index values
# equivalent to `join(ts1, ts2, JoinAll)` call
julia> join(ts1, ts2)
(30 x 2) TS with Date Index

 Index       x1               x1_1
 Date        Float64?         Float64?
────────────────────────────────────────
 2017-01-01        0.768448   0.768448
 2017-01-02        0.940515   0.940515
 2017-01-03        0.673959   0.673959
 2017-01-04        0.395453   0.395453
 2017-01-05        0.313244   0.313244
 2017-01-06        0.662555   0.662555
 2017-01-07        0.586022   0.586022
 2017-01-08        0.0521332  0.0521332
 2017-01-09        0.26864    0.26864
 2017-01-10        0.108871   0.108871
     ⋮              ⋮             ⋮
 2017-01-22  missing          0.291978
 2017-01-23  missing          0.281066
 2017-01-24  missing          0.792931
 2017-01-25  missing          0.20923
 2017-01-26  missing          0.918165
 2017-01-27  missing          0.614255
 2017-01-28  missing          0.802665
 2017-01-29  missing          0.555668
 2017-01-30  missing          0.940782
                         11 rows omitted

# alias to `join()`
julia> cbind(ts1, ts2);

# join only the common index values
julia> join(ts1, ts2, JoinBoth)
(10 x 2) TS with Date Index

 Index       x1         x1_1
 Date        Float64    Float64
──────────────────────────────────
 2017-01-01  0.768448   0.768448
 2017-01-02  0.940515   0.940515
 2017-01-03  0.673959   0.673959
 2017-01-04  0.395453   0.395453
 2017-01-05  0.313244   0.313244
 2017-01-06  0.662555   0.662555
 2017-01-07  0.586022   0.586022
 2017-01-08  0.0521332  0.0521332
 2017-01-09  0.26864    0.26864
 2017-01-10  0.108871   0.108871

# keep index values of `ts1`
julia> join(ts1, ts2, JoinLeft)
(10 x 2) TS with Date Index

 Index       x1         x1_1
 Date        Float64    Float64?
──────────────────────────────────
 2017-01-01  0.768448   0.768448
 2017-01-02  0.940515   0.940515
 2017-01-03  0.673959   0.673959
 2017-01-04  0.395453   0.395453
 2017-01-05  0.313244   0.313244
 2017-01-06  0.662555   0.662555
 2017-01-07  0.586022   0.586022
 2017-01-08  0.0521332  0.0521332
 2017-01-09  0.26864    0.26864
 2017-01-10  0.108871   0.108871

# keep index values of `ts2`
julia> join(ts1, ts2, JoinRight)
(30 x 2) TS with Date Index

 Index       x1               x1_1
 Date        Float64?         Float64
────────────────────────────────────────
 2017-01-01        0.768448   0.768448
 2017-01-02        0.940515   0.940515
 2017-01-03        0.673959   0.673959
 2017-01-04        0.395453   0.395453
 2017-01-05        0.313244   0.313244
 2017-01-06        0.662555   0.662555
 2017-01-07        0.586022   0.586022
 2017-01-08        0.0521332  0.0521332
 2017-01-09        0.26864    0.26864
 2017-01-10        0.108871   0.108871
     ⋮              ⋮             ⋮
 2017-01-22  missing          0.291978
 2017-01-23  missing          0.281066
 2017-01-24  missing          0.792931
 2017-01-25  missing          0.20923
 2017-01-26  missing          0.918165
 2017-01-27  missing          0.614255
 2017-01-28  missing          0.802665
 2017-01-29  missing          0.555668
 2017-01-30  missing          0.940782
                         11 rows omitted

```
"""
function Base.join(ts1::TS, ts2::TS)
    join(ts1, ts2, JoinAll)
end

function Base.join(ts1::TS, ts2::TS, ::Type{JoinBoth})
    result = DataFrames.innerjoin(ts1.coredata, ts2.coredata, on=:Index, makeunique=true)
    return TS(result)
end
Base.join(ts1::TS, ts2::TS, ::Type{JoinInner}) = Base.join(ts1, ts2, JoinBoth)

function Base.join(ts1::TS, ts2::TS, ::Type{JoinAll})
    result = DataFrames.outerjoin(ts1.coredata, ts2.coredata, on=:Index, makeunique=true)
    return TS(result)
end
Base.join(ts1::TS, ts2::TS, ::Type{JoinOuter}) = Base.join(ts1, ts2, JoinAll)

function Base.join(ts1::TS, ts2::TS, ::Type{JoinLeft})
    result = DataFrames.leftjoin(ts1.coredata, ts2.coredata, on=:Index, makeunique=true)
    return TS(result)
end

function Base.join(ts1::TS, ts2::TS, ::Type{JoinRight})
    result = DataFrames.rightjoin(ts1.coredata, ts2.coredata, on=:Index, makeunique=true)
    return TS(result)
end
# alias
cbind = join


text:
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#
# It is included in file Modia3D/Composition/sceneProperties.jl
# in order to be used as default for contact detection in SceneProperties(..)
#
export ContactPair
export NoContactPair
export updateContactPair!
export updateNoContactPair!

using DataStructures

const PairID = Int64

"""
    key = DistanceKey(contact, distanceWithHysteresis, pairID)

Generate a new distance key.
"""
struct DistanceKey
    contact::Bool
    distanceWithHysteresis::Float64
    pairID::PairID
end

"""
    pair = ContactPair(contactPoint1,contactPoint2,contactNormal,obj1,obj2,
                       distanceWithHysteresis)

Generate a new `ContactPair` object of two objects that have contact=true.
"""
mutable struct ContactPair
    contactPoint1::SVector{3,Float64}
    contactPoint2::SVector{3,Float64}
    contactNormal::SVector{3,Float64}
    obj1::Object3D
    obj2::Object3D
    distanceWithHysteresis::Float64

    contactPairMaterial::Modia3D.AbstractContactPairMaterial  # only if contact = true, otherwise not defined

    ContactPair(contactPoint1::SVector{3,Float64}, contactPoint2::SVector{3,Float64}, contactNormal::SVector{3,Float64},
                obj1::Object3D, obj2::Object3D, distanceWithHysteresis::Float64) =
        new(contactPoint1, contactPoint2, contactNormal, obj1, obj2, distanceWithHysteresis)
end

"""
    pair = NoContactPair(contactPoint1,contactPoint2,contactNormal,obj1,obj2,
                         distanceWithHysteresis)

Generate a new `NoContactPair` object of two objects that have contact=false.
"""
mutable struct NoContactPair
    contactPoint1::Union{SVector{3,Float64},Nothing}
    contactPoint2::Union{SVector{3,Float64},Nothing}
    contactNormal::Union{SVector{3,Float64},Nothing}
    obj1::Object3D
    obj2::Object3D
    distanceWithHysteresis::Float64

    contactPairMaterial::Modia3D.AbstractContactPairMaterial  # only if contact = true, otherwise not defined

    NoContactPair(contactPoint1, contactPoint2, contactNormal,
                  obj1::Object3D, obj2::Object3D, distanceWithHysteresis::Float64) =
        new(contactPoint1, contactPoint2, contactNormal, obj1, obj2, distanceWithHysteresis)
end

function updateContactPair!(pair::ContactPair, contactPoint1::SVector{3,Float64}, contactPoint2::SVector{3,Float64}, contactNormal::SVector{3,Float64},
                            obj1::Object3D, obj2::Object3D, distanceWithHysteresis::Float64)::Nothing
    pair.contactPoint1 = contactPoint1
    pair.contactPoint2 = contactPoint2
    pair.contactNormal = contactNormal
    pair.obj1          = obj1
    pair.obj2          = obj2
    pair.distanceWithHysteresis = distanceWithHysteresis
    return nothing
end

function updateNoContactPair!(pair::NoContactPair, contactPoint1, contactPoint2, contactNormal,
                              obj1::Object3D, obj2::Object3D, distanceWithHysteresis::Float64)::Nothing
    pair.contactPoint1 = contactPoint1
    pair.contactPoint2 = contactPoint2
    pair.contactNormal = contactNormal
    pair.obj1          = obj1
    pair.obj2          = obj2
    pair.distanceWithHysteresis = distanceWithHysteresis
    return nothing
end

Base.:isless( value1::ContactPair, value2::ContactPair) = value1.distanceWithHysteresis < value2.distanceWithHysteresis
Base.:isequal(value1::ContactPair, value2::ContactPair) = value1.distanceWithHysteresis == value2.distanceWithHysteresis
Base.:isless( value1::NoContactPair, value2::NoContactPair) = value1.distanceWithHysteresis < value2.distanceWithHysteresis
Base.:isequal(value1::NoContactPair, value2::NoContactPair) = value1.distanceWithHysteresis == value2.distanceWithHysteresis

Base.:isless(key1::DistanceKey, key2::DistanceKey) =   key1.contact &&  key2.contact ? key1.pairID < key2.pairID :
                                                     ( key1.contact && !key2.contact ? true                        :
                                                     (!key1.contact &&  key2.contact ? false                       :
                                                     ( key1.distanceWithHysteresis != key2.distanceWithHysteresis   ?
                                                       key1.distanceWithHysteresis < key2.distanceWithHysteresis    :
                                                       key1.pairID                < key2.pairID)))

Base.:isequal(key1::DistanceKey, key2::DistanceKey) = key1.pairID == key2.pairID

"""
    handler = ContactDetectionMPR_handler(;tol_rel = 1e-4, niter_max=100, neps=sqrt(eps()))

Generate a new contact handler for usage of the MPR algorithm
The handler instance contains all information
about the contact situation.

# Arguments

- `tol_rel`: Relative tolerance to compute the contact point (> 0.0)
- `niter_max`: Maximum number of iterations of the MPR algorithm. If this number is reached,
               an error occurs (> 0).
- `neps`: Small number used to check whether a floating number is close to zero (> 0.0).

"""
mutable struct ContactDetectionMPR_handler <: Modia3D.AbstractContactDetection
  contactPairs::Composition.ContactPairs
  distanceComputed::Bool

  lastContactDict::Dict{PairID,ContactPair}
  contactDict::Dict{    PairID,ContactPair}
  noContactMinVal::Float64

  tol_rel::Float64
  niter_max::Int
  neps::Float64

  # Visualization options
  visualizeContactPoints::Bool
  visualizeSupportPoints::Bool
  defaultContactSphereDiameter::Float64

  function ContactDetectionMPR_handler(;tol_rel   = 1e-4,
                                        niter_max = 100 ,
                                        neps      = sqrt(eps()))
    @assert(tol_rel > 0.0)
    @assert(niter_max > 0)
    @assert(neps > 0.0)
    handler = new()
    handler.distanceComputed = false
    handler.lastContactDict  = Dict{PairID,ContactPair}()
    handler.contactDict      = Dict{PairID,ContactPair}()
    handler.noContactMinVal  = 42.0
    handler.tol_rel          = tol_rel
    handler.niter_max        = niter_max
    handler.neps             = neps
    return handler
  end
end


text:
function simulate_once!(agt::Vector, bt::Billiard{T}, Δt::T) where {T<:AbstractFloat}
    broadcast = ()
    mint = Δt
    minw = 0
    k = 0
    for i = 1:length(agt)
        toc, wall = next_collision(agt[i].p, bt)
        if toc ≤ mint
            mint = toc
            minw = wall
            k = i
        end
    end

    if minw > 0 && typeof(bt[minw]) <: VirtualGate && agt[k].u[bt[minw].color] == 0
        propagate!(agt[k].p, mint + sqrt(eps(mint)))
        broadcast = detected_crossing!(agt[k], bt[minw].color)
    elseif minw > 0
        evolve!(agt[k].p, bt, 1)
    end

    for i = 1:length(agt)
        if i != k
            propagate!(agt[i].p, mint)
            if !isempty(broadcast)
                received_broadcast!(agt[i], broadcast)
            end
        end
    end
end

max_stars_repo_name: gbarsih/Representations-in-Robotics

max_issues_repo_name: gbarsih/Representations-in-Robotics

max_forks_repo_name: gbarsih/Representations-in-Robotics