text:
using ContinuousPOMDPTreeSearchExperiments
using ParticleFilters
using ARDESPOT
using BasicPOMCP
using POMCPOW
using POMDPs
using DiscreteValueIteration
using QMDP
using MCTS
using SubHunt
using POMDPToolbox
using LaserTag
using DataFrames
using CSV
using ProgressMeter
using PmapProgressMeter

@everywhere using LaserTag
@everywhere using ContinuousPOMDPTreeSearchExperiments
@everywhere using POMDPs
@everywhere using ParticleFilters

file_contents = readstring(@__FILE__())

@show max_time = 1.0
@show max_depth = 90
@show exploration = 26.0

P = typeof(gen_lasertag(rng=MersenneTwister(18), robot_position_known=false))

solvers = Dict{String, Union{Solver,Policy}}(
    "qmdp" => QMDPSolver(),

    #=
    Mon 04 Dec 2017 05:11:41 PM PST

    mean(combined[:mean_reward]) = -9.947048511007113
    iteration 88
    mean(d) = [26.3468, 3.62511, 35.0585]
    ev = eigvals(cov(d)) = [2.93607e-11, 6.44598e-8, 0.00381188]
    (eigvecs(cov(d)))[:, j] = [0.00247583, -0.998464, -0.0553574]
    (eigvecs(cov(d)))[:, j] = [0.381521, -0.0502272, 0.922995]
    (eigvecs(cov(d)))[:, j] = [-0.924357, -0.0234052, 0.38081]
    =#

    "pomcpow" => begin
        rng = MersenneTwister(13)
        solver = POMCPOWSolver(tree_queries=10_000_000,
                               criterion=MaxUCB(exploration),
                               final_criterion=MaxTries(),
                               max_depth=max_depth,
                               max_time=max_time,
                               enable_action_pw=false,
                               k_observation=4.0,
                               alpha_observation=1/35,
                               estimate_value=FOValue(ValueIterationSolver()),
                               check_repeat_obs=false,
                               tree_in_info=false,
                               default_action=LaserTag.TAG_ACTION,
                               rng=rng
                              )
    end,

    "pft" => begin
        rng = MersenneTwister(13)
        m = 20
        solver = DPWSolver(n_iterations=typemax(Int),
                           exploration_constant=exploration,
                           depth=max_depth,
                           max_time=max_time,
                           k_state=4.0,
                           alpha_state=1/35,
                           check_repeat_state=false,
                           check_repeat_action=false,
                           tree_in_info=false,
                           estimate_value=RolloutEstimator(QMDPSolver()),
                           enable_action_pw=false,
                           # default_action=ReportWhenUsed(NoGapTag()),
                           rng=rng
                          )
        GBMDPSolver(solver, pomdp->ObsAdaptiveParticleFilter(pomdp,
                                                             LowVarianceResampler(m),
                                                             0.1, rng))
    end,

    "despot" => begin
        rng = MersenneTwister(13)
        # b = IndependentBounds(DefaultPolicyLB(QMDPSolver()), 100.0, check_terminal=true)
        bounds = LaserBounds{P}()
        K = 500
        DESPOTSolver(lambda=0.01,
                     K=K,
                     D=max_depth,
                     max_trials=1_000_000,
                     T_max=max_time,
                     bounds=bounds,
                     bounds_warnings=false,
                     default_action=LaserTag.TAG_ACTION,
                     random_source=MemorizingSource(K, max_depth, rng, min_reserve=10),
                     rng=rng)
    end,

    "d_pomcp" => begin
        rng = MersenneTwister(13)
        POMCPSolver(max_depth=max_depth,
                    max_time=max_time,
                    c=exploration,
                    tree_queries=typemax(Int),
                    tree_in_info=false,
                    default_action=LaserTag.TAG_ACTION,
                    estimate_value=FOValue(ValueIterationSolver()),
                    rng=rng
                   )
    end,

    "pomcpdpw" => begin
        rng = MersenneTwister(13)
        solver = PDPWSolver(tree_queries=10_000_000,
                            c=exploration,
                            max_depth=max_depth,
                            max_time=max_time,
                            enable_action_pw=false,
                            k_observation=4.0,
                            alpha_observation=1/35.0,
                            estimate_value=FOValue(ValueIterationSolver()),
                            check_repeat_obs=false,
                            # default_action=ReportWhenUsed(-1),
                            rng=rng
                           )
    end,
)

@show N=1000

alldata = DataFrame()
for (k, solver) in solvers
# test = ["pft"]
# for (k, solver) in [(s, solvers[s]) for s in test]
# k = "qmdp"
# solver = QMDPSolver()
    @show k
    prog = Progress(N, desc="Creating Simulations...")
    sims = pmap(prog, 1:N) do i
        pomdp = gen_lasertag(rng=MersenneTwister(i+300_000))
        planner = solve(solver, pomdp)
        srand(planner, i+50_000)
        up_rng = MersenneTwister(i+140_000)
        filter = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(100_000), 0.05, up_rng)

        md = Dict(:solver=>k, :i=>i)
        return Sim(pomdp,
            planner,
            filter,
            rng=MersenneTwister(i+70_000),
            max_steps=100,
            metadata=md
           )
    end

    data = run_parallel(sims)
    # data = run(sims)

    rs = data[:reward]
    println(@sprintf("reward: %6.3f ± %6.3f", mean(rs), std(rs)/sqrt(length(rs))))
    if isempty(alldata)
        alldata = data
    else
        alldata = vcat(alldata, data)
    end
end

datestring = Dates.format(now(), "E_d_u_HH_MM")
copyname = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "lasertag_table_$(datestring).jl")
write(copyname, file_contents)
filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "lasertag_$(datestring).csv")
println("saving to $filename...")
CSV.write(filename, alldata)
println("done.")


text:
using LinearAlgebra
using Random
using Printf
using Test

include(joinpath(@__DIR__, "../src/DemandResponse.jl"))
const DR = DemandResponse

import MathOptInterface
const MOI = MathOptInterface

import GLPK
import Linda

# write your own tests here
const testdir = dirname(@__FILE__)

const test_files = [
    # include test file name here (without .jl extension)
    "DER/DER",
    "colgen"
]

for f in test_files
    tp = joinpath(testdir, "$(f).jl")
    include(tp)
end

ext: jl


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -12.548880784274985 <= q <= 15.06837792907793)

                     add_NL_constraint(m, :(gelu(0.4386422028907555 + 0.03049853723546203*gelu(-0.33806544253006265 + -0.5419050015026476*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.7963936693468932*gelu(0.35455126051631725 + -0.3141650319667155*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.19871150900740808*gelu(-0.9712553762218641 + -0.9788234903921222*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) + gelu(0.039843038969352484 + 0.2885956651060697*gelu(-0.33806544253006265 + -0.5419050015026476*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + 0.8477392827604406*gelu(0.35455126051631725 + -0.3141650319667155*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.38122333260621444*gelu(-0.9712553762218641 + -0.9788234903921222*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) + gelu(0.5234531332153054 + -0.26529753551452195*gelu(-0.33806544253006265 + -0.5419050015026476*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.12420368558845318*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.008939690927983968*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + 0.29155449146294243*gelu(0.35455126051631725 + -0.3141650319667155*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + 0.26311214619974344*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + -0.4178451532515153*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))) + -0.357141816330512*gelu(-0.9712553762218641 + -0.9788234903921222*gelu(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])) + -0.6598676192535482*gelu(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])) + 0.7900118209668339*gelu(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --project=@. --color=yes --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

# Not directly used in this script, but needed to enable the glue code in
# GuillotineModels and make it able to use the solver with Requires.jl.
import Gurobi
import TimerOutputs
import CPLEX

import Dates
import Dates: @dateformat_str

# If you want to just check if the calls are being assembled right
# comment the import below and uncomment the mock GuillotineModels
# module after it.
import GuillotineModels
# Mock GuillotineModels.run used to check if the script is correct.
#=
module GuillotineModels
	import TimerOutputs
	const TIMER = TimerOutputs.TimerOutput()

	struct TimeoutError <: Exception end
	function Base.showerror(io :: IO, e :: TimeoutError)
		print(io, "mock timeout error")
	end

	module CommandLine
		import ..TimeoutError
		function run(args; kwargs...)
			@show args
			@show kwargs
			if any(s -> !isnothing(match(r".*/Hchl6s", s)), args)
				throw(TimeoutError())
			end
			if any(s -> !isnothing(match(r".*/A5", s)), args)
				error("mock ErrorException")
			end
		end
	end
end
=#

function save_output_in_path(f, path :: AbstractString)
	open(path, "a+") do file
		save_output_in_file(f, file)
	end
end
function save_output_in_file(f, file :: IO)
  redirect_stdout(file) do
		redirect_stderr(file) do
			f()
			# See: discourse.julialang.org/t/cannot-capture-log-with-jump/55365
			# My thanks to Miles Lubin
			Base.Libc.flush_cstdio()
		end
	end
end

macro display_error(io, e)
	return quote
		println("run_ended_by_exception = true")
		showerror($(esc(io)), $(esc(e)))
		println($(esc(io)))
		show($(esc(io)), "text/plain", stacktrace(catch_backtrace()))
		println($(esc(io)))
	end
end

function safe_run(args, supported_solvers, implemented_models, verbose = true)
	try
		GuillotineModels.CommandLine.run(
			args; supported_solvers = supported_solvers,
			implemented_models = implemented_models
		)
		verbose && println("run_ended_by_exception = false")
	catch e
		# If `verbose == false` this is probably a mock run, and we do want to
		# have exception logging in a mock run because no exceptions should
		# happen inside a mock run (not even timeouts).
		if isa(e, GuillotineModels.TimeoutError)
			@display_error stderr e
		elseif isa(e, LoadError) && isa(e.error, GuillotineModels.TimeoutError)
			@display_error stderr e.error
		else
			@display_error stderr e
			open("UNEXPECTED_EXCEPTION_README.log", "a") do io
				println(io, Dates.format(
					Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
				))
				println(io, args)
				@display_error io e
			end
		end
	end
	# The GC should happen in a mock run, but not the printing.
	extra_gc_time = @elapsed GC.gc()
	if verbose
		@show extra_gc_time
		TimerOutputs.print_timer(GuillotineModels.TIMER; allocations = false)
		println()
	end
	TimerOutputs.reset_timer!(GuillotineModels.TIMER)
	return
end

function saved_run(
	args               :: Vector{String},
	supported_solvers  :: Vector{Symbol},
	implemented_models :: Vector{Symbol},
	output_folder      :: String,
	verbose            :: Bool = true
)
	format = dateformat"yyyy-mm-ddTHH:MM:SS"
	verbose && println("Next run arguments: $args")
	verbose && println("Started at: ", Dates.format(Dates.now(), format))
	filepath, io = mktemp(output_folder; cleanup = false)
	save_output_in_file(io) do
		_, t, bytes, gctime, _ = @timed safe_run(
			args, supported_solvers, implemented_models, true
		)
		println("run_total_time = $t")
		println("run_bytes_allocated = $bytes")
		println("run_time_spent_on_gc = $gctime")
		println("this_data_file = $filepath")
	end
	close(io)
	verbose && println("Finished at: ", Dates.format(Dates.now(), format))
	println("Output saved to: $(filepath)")
	return filepath
end

# The concept of a batch is a set of runs that share the same model, solver,
# and most options but may vary the instances and the random seeds.
# The main utility of a batch is that we solve a mock instance before the
# other instances and then we have the guarantee that everything that depends
# on model, solver, and most options will be already compiled by JIT before
# we make serious time measurements (theoretically, the change of seed and
# instance should not change the inner methods called).
function run_batch(
	problem           :: String
	, instance_format :: String
	, formulation     :: String
	, solver          :: String
	, mock_inst_path  :: String
	, instance_paths  :: AbstractVector{String}
	; options         :: Vector{String} = String[]
	, solver_seeds    :: Vector{Int} = Int[]
	, output_folder   :: String = "./"
)
	supported_solvers = [Symbol(solver)]
	implemented_models = [Symbol(formulation)]
	# The mock should not print anything, and even if it prints something,
	# it will not go into the file (it is not entirely supressed to make
	# the failure of these assumptions visible so we can correct them).
	mock_output_flags = [
		"--no-csv-output", "--$solver-no-output", "--$formulation-quiet"
	]
	tinkered_options = append!(["--warm-jit", "no"], options)
	mock_flags = vcat(tinkered_options, mock_output_flags)
	mock_args = append!(
		[problem, instance_format, formulation, solver, mock_inst_path],
		mock_flags
	)
	safe_run(mock_args, supported_solvers, implemented_models, false)
	for inst_path in instance_paths
		non_seed_args = append!(
			[problem, instance_format, formulation, solver, inst_path],
			tinkered_options
		)
		if isempty(solver_seeds)
			filename = saved_run(
				non_seed_args, supported_solvers, implemented_models, output_folder
			)
		else
			for seed in solver_seeds
				args_with_seed = append!(["--$solver-seed", "$seed"], non_seed_args)
				filename = saved_run(
					args_with_seed, supported_solvers, implemented_models, output_folder
				)
			end
		end
	end
	return
end

#=
function run_experiments(
	solver = "CPLEX"
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	all_instance_names = vcat(
		"gcut" .* string.(collect(2:12))
		, ["STS4", "STS4s", "A5"]
		, "Hchl" .* split("2 3s 4s 6s 7s")
		, "CW" .* split("1 2 3")
		, "CHL" .* split("1 1s 6 7")
		, "CU" .* split("1 2")
		, "okp" .* split("1 2 3 4 5")
		, "P1_100_200_25_" .* split("1 2 3 4 5")
	)
	all_instance_paths = instance_folder .* all_instance_names
	easy_instance_names = vcat(
		"okp" .* split("1 4 5")
		, ["CU1", "STS4", "STS4s", "gcut9"]
	)
	easy_instance_paths = instance_folder .* easy_instance_names

	tl = 60.0

	#=
	let # just to avoid the variable names from leaking
		common_flags = ["--do-not-solve", "--PPG2KP-faithful2furini2016"]
		option_sets = [
			String[]
			, ["--PPG2KP-no-redundant-cut", "--PPG2KP-no-cut-position"]
			, ["--PPG2KP-no-redundant-cut"]
			, ["--PPG2KP-no-cut-position"]
			, ["--PPG2KP-round2disc"]
		]
		for options in option_sets
			append!(options, common_flags) # note: this changes `option_sets`
			run_batch(
				"PPG2KP", solver, mock_inst_path, all_instance_paths, options;
				output_path = output_path
			)
		end
	end
	=#

	let # just to avoid the variable names from leaking
		common_options = [
			"--generic-time-limit", "$tl", "--PPG2KP-building-time-limit", "$tl",
			"--PPG2KP-verbose"
		]
		#=option_sets = [
			String[], ["--PPG2KP-faithful2furini2016"], #=["--PPG2KP-no-pricing"],
			["--PPG2KP-no-pricing", "--PPG2KP-faithful2furini2016"]=#
			["--PPG2KP-round2disc"], ["--PPG2KP-round2disc",
			"--PPG2KP-faithful2furini2016"]
			#["--PPG2KP-no-cut-position", "--PPG2KP-no-redundant-cut"]
		]=#
		option_sets = [
			["--PPG2KP-pricing", "none"], # revised, no pricing
			String[], # revised, becker pricing
			["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "none"], # furini, no pricing
			["--PPG2KP-faithful2furini2016"] # furini, furini pricing
		]
		#solver_seeds = [1] #collect(1:10)
		for options in option_sets
			append!(options, common_options) # NOTE: changes `option_sets` elements
			run_batch(
				"PPG2KP", solver, easy_instance_paths[1], easy_instance_paths;
				options = options, #solver_seeds = solver_seeds,
				output_folder = output_folder
			)
		end
	end
end
=#

# The 7 instances that Furini 2016 method without pricing is able to solve
# within their time limit.
const EASY_SEVEN = vcat(
	"okp" .* split("1 4 5"), ["CU1", "STS4", "STS4s", "gcut9"]
)

function run_LP_method_experiment(
	solver = "CPLEX"
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	easy_instance_paths = instance_folder .* EASY_SEVEN

	time_limit = 3600.0

	common_options = [
		"--generic-time-limit", "$time_limit", "--PPG2KP-building-time-limit",
		"$time_limit", "--PPG2KP-verbose", "--PPG2KP-faithful2furini2016"
	]
	option_sets = [
		# No pricing, use dual simplex.
		["--PPG2KP-pricing", "none", "--Gurobi-LP-method", "1"],
		# No pricing, use barrier.
		["--PPG2KP-pricing", "none", "--Gurobi-LP-method", "2"],
		# Furini pricing, use dual simplex.
		["--PPG2KP-pricing", "furini", "--Gurobi-LP-method", "1"],
		# Furini pricing, use barrier.
		["--PPG2KP-pricing", "furini", "--Gurobi-LP-method", "2"],
		# Furini pricing, use barrier outside of the pricing, in the pricing use
		# dual simplex.
		[
			"--PPG2KP-pricing", "furini", "--Gurobi-LP-method", "2",
			"--PPG2KP-Gurobi-LP-method-inside-furini-pricing", "1"
		]
	]
	solver_seeds = [1]
	for options in option_sets
		append!(options, common_options) # NOTE: changes `option_sets` elements
		run_batch(
			# TODO: the mock instance probably should not be hardcoded below
			"PPG2KP", solver, instance_folder * "CU1", easy_instance_paths;
			options = options, solver_seeds = solver_seeds,
			output_folder = output_folder
		)
	end
	return
end

# The 59 instances presented in table 2.1, page 30,
# from DOI: 10.6092/unibo/amsdottorato/7399 (order was kept).
const THOMOPULOS_THESIS_INSTANCES = vcat(
	# unweighted (37)
	"gcut" .* string.(1:12)
	, String.(split("wang20 2s 3s A1s A2s STS2s STS4s"))
	, ["OF1", "OF2", "W", "CHL1s", "CHL2s"]
	, "A" .* string.(3:5)
	, "CHL" .* string.(5:7)
	, ["CU1", "CU2"]
	, "Hchl" .* split("3s 4s 6s 7s 8s")
	# weighted (22)
	, "cgcut" .* string.(1:3)
	, "okp" .* string.(1:5)
	, String.(split("HH 2 3 A1 A2 STS2 STS4 CHL1 CHL2"))
	, "CW" .* string.(1:3)
	, "Hchl" .* ["2", "9"]
) :: Vector{String}

function run_faithful_reimplementation_experiment(
	solver = "CPLEX"
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	#HARD_SIX = ["gcut11", "gcut12", "okp3", "okp2", "Hchl6s", "Hchl7s"]
	instance_paths = instance_folder .* THOMOPULOS_THESIS_INSTANCES
	time_limit = 10800.0 # three hours

	common_options = [
		"--generic-time-limit", "$time_limit", "--PPG2KP-building-time-limit",
		"$time_limit", "--PPG2KP-verbose", "--PPG2KP-faithful2furini2016",
		"--do-not-solve"
	]
	option_sets = [
		# The Complete PP-G2KP (no reductions).
		["--PPG2KP-pricing", "none", "--PPG2KP-no-redundant-cut",
			"--PPG2KP-no-cut-position"],
		# The Complete PP-G2KP (only cut position).
		["--PPG2KP-pricing", "none", "--PPG2KP-no-redundant-cut"],
		# The Complete PP-G2KP (only redundant cut).
		["--PPG2KP-pricing", "none", "--PPG2KP-no-cut-position"],
		# The Complete PP-G2KP (no reductions, no pricing).
		["--PPG2KP-pricing", "none"],
		# The PP-G2KP (i.e., with pricing).
		["--PPG2KP-pricing", "furini",
			"--PPG2KP-Gurobi-LP-method-inside-furini-pricing", "1"]
	]
	solver_seeds = [1]
	for options in option_sets
		append!(options, common_options) # NOTE: changes `option_sets` elements
		# gcut6 was chosen as the mock instance because it has median
		# times in all option sets, so it probably enter the longer
		# computation paths on all of them.
		run_batch(
			"PPG2KP", solver, instance_folder * "gcut6", instance_paths;
			options = options, solver_seeds = solver_seeds,
			output_folder = output_folder
		)
	end
	return
end

function run_comparison_experiment(
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	instance_paths = instance_folder .* THOMOPULOS_THESIS_INSTANCES
	#instance_paths = [instance_folder .* "A5"]
	time_limit = 10800.0 # three hours

	common_options = [
		"--generic-time-limit", "$time_limit", "--PPG2KP-building-time-limit",
		"$time_limit", "--PPG2KP-verbose"
	]
	solver_options = Dict{String, Vector{String}}(
		"CPLEX" => ["--CPLEX-root-relax-method", "barrier"],
		"Gurobi" => [
			"--Gurobi-LP-method", "2" #= 2 == barrier =#,
			"--PPG2KP-Gurobi-LP-method-inside-furini-pricing", "1" #= 1 == dual =#
		]
	)
	option_sets = [
		# Just the revised model.
		#["--PPG2KP-pricing", "none"],
		# The revised model with our reduction.
		#["--PPG2KP-pricing", "none", "--PPG2KP-round2disc"],
		# The revised model with our reduction and warm-start.
		#["--PPG2KP-pricing", "none", "--PPG2KP-round2disc",
		#	"--PPG2KP-MIP-start", "guaranteed"],
		# The revised model with our reduction, pricing, and warm-start.
		#["--PPG2KP-pricing", "furini", "--PPG2KP-round2disc",
		#	"--PPG2KP-MIP-start", "guaranteed"],
		# The revised model with our reduction, pricing, and warm-start,
		# but removal of unreachable disabled.
		#["--PPG2KP-pricing", "furini", "--PPG2KP-round2disc",
		#	"--PPG2KP-MIP-start", "guaranteed",
		#	"--PPG2KP-do-not-purge-unreachable"],
		# The original model (no pricing, just their reductions).
		["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "none"],
		# The original model with our reduction too.
		["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "none",
			"--PPG2KP-round2disc"],
		# The orifinal model with our reduction and warm-start.
		["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "none",
			"--PPG2KP-round2disc", "--PPG2KP-MIP-start", "guaranteed"],
		# The original model with furini pricing (Priced PPG2KP).
		#["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "furini"]
		# The original model with our reduction plus furini pricing.
		#["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "furini",
		#	"--PPG2KP-round2disc"]
	]
	solver_seeds = [1]#, 2, 3]
	for solver in [#="CPLEX",=# "Gurobi"]
		for options in option_sets
			append!(options, common_options) # NOTE: changes `option_sets` elements
			# gcut6 was chosen as the mock instance because it has median
			# times in all option sets, so it probably enter the longer
			# computation paths on all of them.
			run_batch(
				"PPG2KP", solver,
				instance_folder * "gcut6",
				instance_paths;
				options = vcat(options, solver_options[solver]),
				solver_seeds = solver_seeds,
				output_folder = output_folder
			)
		end
	end
	return
end

const VELASCO_AND_UCHOA_80 = String[
	"P$(i)_$(L)_$(W)_$(n)_$(s)" for i in 1:4 for (L, W) in
	([(100, 200), (100, 400)], [(200, 100), (400, 100)],
		[(150, 150), (250, 250)], [(150, 150), (250, 250)])[i] for n in (25, 50)
	for s in 1:5
]

function run_vel_uchoa_experiment(
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	instance_paths = instance_folder .* VELASCO_AND_UCHOA_80
	time_limit = 10800.0 # three hours

	common_options = [
		"--generic-time-limit", "$time_limit", "--PPG2KP-building-time-limit",
		"$time_limit", "--PPG2KP-verbose"
	]
	solver_options = Dict{String, Vector{String}}(
		#"CPLEX" => ["--CPLEX-root-relax-method", "barrier"],
		"Gurobi" => [
			"--Gurobi-LP-method", "2" #= 2 == barrier =#,
			"--PPG2KP-Gurobi-LP-method-inside-furini-pricing", "1" #= 1 == dual =#
		]
	)
	option_sets = [
		# Just the revised model.
		#["--PPG2KP-pricing", "none"],
		# The revised model with our reduction.
		#["--PPG2KP-pricing", "none", "--PPG2KP-round2disc"],
		# The revised model with our reduction, pricing, and warm-start.
		["--PPG2KP-pricing", "furini", "--PPG2KP-round2disc",
			"--PPG2KP-MIP-start", "guaranteed"],
		# The revised model with our reduction and warm-start.
		["--PPG2KP-pricing", "none", "--PPG2KP-round2disc",
			"--PPG2KP-MIP-start", "guaranteed"],
		# The revised model with our reduction, pricing, and warm-start,
		# but removal of unreachable disabled.
		#["--PPG2KP-pricing", "furini", "--PPG2KP-round2disc",
		#	"--PPG2KP-MIP-start", "guaranteed",
		#	"--PPG2KP-do-not-purge-unreachable"],
		# The original model (no pricing, just their reductions).
		#["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "none"],
		# The original model with our reduction too.
		#["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "none",
		#	"--PPG2KP-round2disc"],
		# The orifinal model with our reduction and warm-start.
		#["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "none",
		#	"--PPG2KP-round2disc", "--PPG2KP-MIP-start", "guaranteed"],
		# The original model with furini pricing (Priced PPG2KP).
		#["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "furini"]
		# The original model with our reduction plus furini pricing.
		#["--PPG2KP-faithful2furini2016", "--PPG2KP-pricing", "furini",
		#	"--PPG2KP-round2disc"]
	]
	solver_seeds = [1]#, 2, 3]
	for solver in [#="CPLEX",=# "Gurobi"]
		for options in option_sets
			append!(options, common_options) # NOTE: changes `option_sets` elements
			# gcut6 was chosen as the mock instance because we have no idea yet which
			# instance from the velasco and uchoa is a good mock
			run_batch(
				"PPG2KP", solver,
				instance_folder * "gcut6",
				instance_paths;
				options = vcat(options, solver_options[solver]),
				solver_seeds = solver_seeds,
				output_folder = output_folder
			)
		end
	end
	return
end

# DATASET C OF https://doi.org/10.1016/j.ejor.2010.01.039
const DATASET_C = vcat(
	"A" .* string.(1:5),
	"CHL" .* string.([1, 2, 5, 6, 7]),
	"CU" .* ["1", "2"],
	"CW" .* string.(1:3),
	"Hchl" .* ["2", "9"],
	"Hchl" .* string.([3, 4, 6, 7, 8]) .* "s",
	"HH",
	"OF" .* ["1", "2"],
	"STS" .* ["2", "4"],
	"W",
	"2",
	"3"
)

# DATASET D OF https://doi.org/10.1016/j.ejor.2010.01.039
const DATASET_D = "ATP" .* string.(30:49)

const GCUTS = "gcut" .* string.(1:12)

function run_lagos_experiment(
	all_instances :: AbstractVector{String}
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	instance_paths = instance_folder .* all_instances
	time_limit = 10800.0 # three hours

	common_options = [
		"--generic-time-limit", "$time_limit", "--PPG2KP-building-time-limit",
		"$time_limit", "--PPG2KP-verbose"
	]
	option_sets = [
		# The revised model, multiple-threads.
		[ # FOR YOJIMBO
			"--PPG2KP-pricing", "none", "--PPG2KP-round2disc",
			"--Gurobi-threads", "12", "--Gurobi-LP-method", "3", # 3 == concurrent
    ],
		# The revised model, single-thread.
		[ # FOR YOJIMBO
			"--PPG2KP-pricing", "none", "--PPG2KP-round2disc",
			"--Gurobi-LP-method", "2", # 2 == barrier
		],
		# The original model, multiple-threads.
		[ # FOR LEVIATHAN
			"--PPG2KP-pricing", "none",
			"--PPG2KP-faithful2furini2016",
			"--Gurobi-threads", "12", "--Gurobi-LP-method", "3", # 3 == concurrent
    ],
		# The original model, single-thread.
		[ # FOR LEVIATHAN
			"--PPG2KP-pricing", "none",
			"--PPG2KP-faithful2furini2016",
			"--Gurobi-LP-method", "2", # 2 == barrier
		],
		# The original model with plate-size normalization, multiple-threads.
		[ # FOR ODIN
			"--PPG2KP-pricing", "none", "--PPG2KP-round2disc",
			"--PPG2KP-faithful2furini2016",
			"--Gurobi-threads", "12", "--Gurobi-LP-method", "3", # 3 == concurrent
    ],
		# The original model with plate-size normalization, single-thread.
		[ # FOR ODIN
			"--PPG2KP-pricing", "none", "--PPG2KP-round2disc",
			"--PPG2KP-faithful2furini2016",
			"--Gurobi-LP-method", "2", # 2 == barrier
		],
	]
	solver_seeds = [1]
	for options in option_sets
		append!(options, common_options) # NOTE: changes `option_sets` elements
		run_batch(
			"PPG2KP", "Gurobi",
			instance_folder * "STS4", # This is the mock instance.
			instance_paths;
			options = options,
			solver_seeds = solver_seeds,
			output_folder = output_folder
		)
	end
	return
end

const CWs = "CW" .* string.(1:11)
const CUs = "CU" .* string.(1:11)

function run_rotation_experiment(
	solver :: String,
	all_instances :: AbstractVector{String},
	mock_instance :: String,
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	instance_paths = instance_folder .* all_instances
	time_limit = 7200.0 # two hours

	common_options = [
		"--generic-time-limit", "$time_limit", "--PPG2KP-building-time-limit",
		"$time_limit", "--PPG2KP-verbose", "--PPG2KP-round2disc",
		"--PPG2KP-pricing", "none"
	]
	option_sets = [
		["--PPG2KP-allow-rotation"]
	]
	@assert solver in ("Gurobi", "CPLEX")
	if solver == "Gurobi"
		append!.(option_sets, (["--Gurobi-LP-method", "2"],)) # barrier
	else
		append!.(option_sets, (["--CPLEX-LP-method", "4"],)) # barrier
	end
	solver_seeds = [1]
	for options in option_sets
		append!(options, common_options) # NOTE: changes `option_sets` elements
		run_batch(
			"G2KP", "Classic_G2KP", "PPG2KP", solver,
			instance_folder * mock_instance, # This is the mock instance.
			instance_paths;
			options = options,
			solver_seeds = solver_seeds,
			output_folder = output_folder
		)
	end
	return
end

function save_models(
	all_instances :: AbstractVector{String},
	mock_instance :: String,
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	),
	extension = "mps",
	mps_folder = joinpath("./$extension/", last(splitpath(instance_folder)))
)
	# extension should be ascii utf8 so this should be safe
	@assert extension[1] != '.'
	isdir(output_folder) || mkpath(output_folder)
	isdir(mps_folder) || mkpath(mps_folder)
	instance_paths = joinpath.((instance_folder,), all_instances)

	common_options = [
		"--PPG2KP-verbose", "--PPG2KP-pricing", "none",
		"--do-not-solve", "--save-model", joinpath(
			mps_folder, "enhanced2-\$<instance_file>" *
			"-\$<PPG2KP-allow-rotation>.$extension"
		)
	]
	option_sets = [
		["--PPG2KP-round2disc", "--PPG2KP-allow-rotation",
			"--PPG2KP-mirror-plate"],
		#["--PPG2KP-faithful2furini2016", "--PPG2KP-allow-rotation",
		#	"--PPG2KP-mirror-plate"],
		["--PPG2KP-round2disc"], # no rotation
		#["--PPG2KP-faithful2furini2016"], # no rotation
	]
	for options in option_sets
		append!(options, common_options) # NOTE: changes `option_sets` elements
		run_batch(
			"G2KP", "Simple_CPG_SLOPP", "PPG2KP", "NoSolver",
			instance_folder * mock_instance, # This is the mock instance.
			instance_paths; options = options, output_folder = output_folder
		)
	end
	return
end

const HARD4 = String[
	"Hchl4s", "Hchl7s", "okp2", "okp3"
]

function run_hybridization_experiment(
	solver :: String,
	all_instances :: AbstractVector{String},
	mock_instance :: String,
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	instance_paths = instance_folder .* all_instances
	time_limit = 3600.0 * 3 # three hours

	common_options = [
		"--generic-time-limit", "$time_limit", "--PPG2KP-building-time-limit",
		"$time_limit", "--PPG2KP-verbose", "--PPG2KP-round2disc",
		"--PPG2KP-pricing", "none"
	]
	option_sets = [
		String[], # Only common options
		["--PPG2KP-hybridize-with-restricted"],
		["--PPG2KP-hybridize-with-restricted",
			"--PPG2KP-aggressive-hybridization"],
	]
	@assert solver in ("Gurobi", "CPLEX")
	if solver == "Gurobi"
		append!.(option_sets, (["--Gurobi-LP-method", "2"],)) # barrier
	else
		append!.(option_sets, (["--CPLEX-LP-method", "4"],)) # barrier
	end
	solver_seeds = [1]
	for options in option_sets
		append!(options, common_options) # NOTE: changes `option_sets` elements
		run_batch(
			"G2KP", "Classic_G2KP", "PPG2KP", solver,
			instance_folder * mock_instance, # This is the mock instance.
			instance_paths;
			options = options,
			solver_seeds = solver_seeds,
			output_folder = output_folder
		)
	end
	return
end

const Clautiaux42 = String[
	"E00N10", "E00N15", "E00N23", "E00X23", "E02F17", "E02F20", "E02F22",
	"E02N20", "E03N10", "E03N15", "E03N16", "E03N17", "E03X18", "E04F15",
	"E04F17", "E04F19", "E04F20", "E04N15", "E04N17", "E04N18", "E05F15",
	"E05F18", "E05F20", "E05N15", "E05N17", "E05X15", "E07F15", "E07N10",
	"E07N15", "E07X15", "E08F15", "E08N15", "E10N10", "E10N15", "E10X15",
	"E13N10", "E13N15", "E13X15", "E15N10", "E15N15", "E20F15", "E20X15"
]

const HopperTurton_C = String[
	"c$(c)-p$(p)" for c in 1:7 for p in 1:3
]

function run_G2OPP_exploratory_experiment(
	solver :: String,
	all_instances :: AbstractVector{String},
	mock_instance :: String,
	; instance_folder :: String = "../instances/"
	, output_folder   :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
)
	isdir(output_folder) || mkpath(output_folder)
	instance_paths = instance_folder .* all_instances
	time_limit = 3600.0 * 3 # three hours

	common_options = [
		"--generic-time-limit", "$time_limit", "--PPG2KP-building-time-limit",
		"$time_limit", "--PPG2KP-verbose", "--PPG2KP-round2disc",
		"--PPG2KP-pricing", "none",	]
	option_sets = [
		String["--Gurobi-raw-parameter",
			"Pair{String, Any}[\"BarHomogeneous\" => 1]"],
		String["--Gurobi-raw-parameter",
			"Pair{String, Any}[\"BarHomogeneous\" => 1, \"Cutoff\" => 1.5]"],
	]
	@assert solver in ("Gurobi", "CPLEX")
	if solver == "Gurobi"
		append!.(option_sets, (["--Gurobi-LP-method", "2"],)) # barrier
	else
		append!.(option_sets, (["--CPLEX-LP-method", "4"],)) # barrier
	end
	solver_seeds = [1]
	for options in option_sets
		append!(options, common_options) # NOTE: changes `option_sets` elements
		run_batch(
			"G2CSP", "CPG_SSSCSP", "PPG2KP", solver,
			instance_folder * mock_instance, # This is the mock instance.
			instance_paths;
			options = options,
			solver_seeds = solver_seeds,
			output_folder = output_folder
		)
	end
	return
end

# The ordering of the instances is the ideal for prioritising solving
# a representative subset of them: Seed > Size > Category.
const CLASS = [
	"cl_$(lpad(c, 2, '0'))_$(lpad(n, 3, '0'))_$(lpad(s, 2, '0'))"
	for s in 1:10 for n in 20:20:100 for c in 1:10
]

const CLASS_50 = [
	"cl_$(lpad(c, 2, '0'))_$(lpad(n, 3, '0'))_01"
	for n in 20:20:100 for c in 1:10
]

const A = [
	"A42",  "A4",  "A8",  "A3",  "A6", "A39", "A43",  "A1",
	"A22",  "A2", "A30", "A10", "A12", "A34", "A36",  "A7",
	"A26", "A31", "A41", "A37", "A28", "A17",  "A5", "A35",
	"A25", "A38", "A40", "A27", "A32", "A20", "A24", "A29",
	"A23", "A13", "A33", "A15", "A21", "A19", "A18", "A11",
	"A16",  "A9", "A14"
]

const A_MKP = [
	"A02_M2", "A02_M4", "A02_M8", "A03_M2", "A05_M2", "A05_M4", "A07_M2",
	"A07_M4", "A09_M2", "A09_M4", "A09_M8", "A11_M2", "A11_M4", "A11_M8",
	"A12_M2", "A12_M4", "A13_M2", "A13_M4", "A14_M2", "A14_M4", "A14_M8",
	"A15_M2", "A15_M4", "A15_M8", "A16_M2", "A16_M4", "A16_M8", "A17_M2",
	"A18_M2", "A18_M4", "A18_M8", "A19_M2", "A19_M4", "A19_M8", "A20_M2",
	"A20_M4", "A20_M8", "A21_M2", "A21_M4", "A21_M8", "A23_M2", "A23_M4",
	"A24_M2", "A24_M4", "A24_M8", "A25_M2", "A25_M4", "A26_M2", "A27_M2",
	"A27_M4", "A27_M8", "A28_M2", "A28_M4", "A29_M2", "A29_M4", "A29_M8",
	"A31_M2", "A32_M2", "A32_M4", "A32_M8", "A33_M2", "A33_M4", "A33_M8",
	"A34_M2", "A35_M2", "A35_M4", "A36_M2", "A37_M2", "A38_M2", "A38_M4",
	"A38_M8", "A40_M2", "A40_M4", "A41_M2", "A41_M4", "A42_M2", "A43_M2"
]

const CW_MKP = [
	"CW01_M2", "CW01_M4", "CW02_M2", "CW02_M4", "CW03_M2", "CW03_M4",
	"CW04_M2", "CW04_M4", "CW05_M2", "CW05_M4", "CW06_M2", "CW06_M4",
	"CW06_M8", "CW07_M2", "CW07_M4", "CW07_M8", "CW08_M2", "CW08_M4",
	"CW08_M8", "CW09_M2", "CW09_M4", "CW09_M8", "CW10_M2", "CW10_M4",
	"CW10_M8", "CW11_M2", "CW11_M4", "CW11_M8"
]

function run(
	solver          :: String,
	problem         :: String,
	format          :: String,
	all_instances   :: AbstractVector{String},
	mock_instance   :: String,
	instance_folder :: String
	; output_folder :: String = "./experiments_outputs/" * Dates.format(
		Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS"
	)
	, mock_instance_folder = instance_folder
	, MIPGap :: Float64 = 1e-6
	, time_limit :: Float64 = 3600.0
	, Gurobi_LP_method = 2
	, CPLEX_LP_method = 4
	, Gurobi_raw_parameters = String[]
	, CPLEX_raw_parameters = String[]
	, solver_seeds = (1,)
	, option_sets :: Vector{Vector{String}} = [["--PPG2KP-round2disc"]]
)
	isdir(output_folder) || mkpath(output_folder)
	instance_paths = instance_folder .* all_instances

	option_sets = deepcopy(option_sets)

	common_options = [
		"--generic-time-limit", "$time_limit",
		"--PPG2KP-building-time-limit", "$time_limit",
		"--PPG2KP-verbose",
		"--PPG2KP-pricing", "none",
	]
	@assert solver in ("Gurobi", "CPLEX")
	if solver == "Gurobi"
		default_gurobi_raw_params = [
			"\"NumericFocus\" => 3",
			"\"MIPGap\" => $(MIPGap)",
		]
		all_gurobi_raw_params = vcat(
			default_gurobi_raw_params, Gurobi_raw_parameters
		)
		gurobi_raw_params_str = join(string.(all_gurobi_raw_params), ',')
		append!.(option_sets, ([
			"--Gurobi-raw-parameters", "Pair{String, Any}[$(gurobi_raw_params_str)]"
		],))
		append!.(option_sets, (["--Gurobi-LP-method", string(Gurobi_LP_method)],))
	else
		append!.(option_sets, ([
			"--CPLEX-raw-parameters",
			"Pair{String, Any}[\"CPXPARAM_MIP_Tolerances_MIPGap\" => $(MIPGap)]"
		],))
		append!.(option_sets, (["--CPLEX-LP-method", string(CPLEX_LP_method)],))
	end

	for options in option_sets
		append!(options, common_options) # NOTE: changes `option_sets` elements
		run_batch(
			problem, format, "PPG2KP", solver,
			joinpath(mock_instance_folder, mock_instance), # This is the mock instance.
			instance_paths;
			options = options,
			solver_seeds = collect(solver_seeds),
			output_folder = output_folder
		)
	end
	return
end

#=
run(
	"Gurobi", "G2CSP", "CPG_SSSCSP", CLASS, CLASS[1],
	"../instances/G2CSP/CLASS/"; MIPGap = 1e-4
)

run(
	"Gurobi", "G2CSP", "CPG_SSSCSP", A, A[1], "../instances/G2CSP/A/";
	MIPGap = 1e-4
)

run(
	"Gurobi", "G2MKP", "CPG_MHLOPPW", A_MKP, A_MKP[1],
	"../instances/G2MKP/"; MIPGap = 1e-8
)

run(
	"Gurobi", "G2MKP", "CPG_MHLOPPW", CW_MKP, CW_MKP[1],
	"../instances/G2MKP/"; MIPGap = 1e-7
)

run(
	"Gurobi", "G2OPP", "CPG_SSSCSP", Clautiaux42, Clautiaux42[1],
	"../instances/G2OPP/Clautiaux42/"
)

run(
	"Gurobi", "G2OPP", "CPG_SSSCSP", HopperTurton_C, HopperTurton_C[1],
	"../instances/G2OPP/HopperTurton/C/"
)
=#

run(
	"Gurobi", "G2CSP", "CPG_SSSCSP", [CLASS[1]], CLASS[1],
	"../instances/G2CSP/CLASS/"; MIPGap = 1e-4
)

run(
	"Gurobi", "G2CSP", "CPG_SSSCSP", [A[1]], A[1], "../instances/G2CSP/A/";
	MIPGap = 1e-4
)

run(
	"Gurobi", "G2MKP", "CPG_MHLOPPW", [A_MKP[1]], A_MKP[1],
	"../instances/G2MKP/"; MIPGap = 1e-8
)

run(
	"Gurobi", "G2MKP", "CPG_MHLOPPW", [CW_MKP[1]], CW_MKP[1],
	"../instances/G2MKP/"; MIPGap = 1e-7
)

run(
	"Gurobi", "G2OPP", "CPG_SSSCSP", [Clautiaux42[1]], Clautiaux42[1],
	"../instances/G2OPP/Clautiaux42/"
)

run(
	"Gurobi", "G2OPP", "CPG_SSSCSP", [HopperTurton_C[1]], HopperTurton_C[1],
	"../instances/G2OPP/HopperTurton/C/"
)


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Special Orthogonal group" begin
    for n in [2, 3]
        G = SpecialOrthogonal(n)
        @test repr(G) == "SpecialOrthogonal($n)"
        M = base_manifold(G)
        @test M === Rotations(n)
        p = Matrix(I, n, n)
        @test is_default_metric(MetricManifold(G, EuclideanMetric()))

        types = [Matrix{Float64}]

        if n == 2
            ω = [[1.0], [-1.0], [3.0]]
            pts = [exp(M, p, hat(M, p, ωi)) for ωi in ω]
            Xpts = [hat(M, p, [-1.5]), hat(M, p, [0.5])]
        else
            ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
            pts = [exp(M, p, hat(M, p, ωi)) for ωi in ω]
            Xpts = [hat(M, p, [-1.0, 2.0, 0.5]), hat(M, p, [1.0, 0.0, 0.5])]
        end

        ge = allocate(pts[1])
        identity_element!(G, ge)
        @test isapprox(ge, I; atol=1e-10)

        gI = Identity(G)
        gT = allocate_result(G, exp, gI, log(G, pts[1], pts[2]))
        @test size(gT) == size(ge)
        @test eltype(gT) == eltype(ge)
        gT = allocate_result(G, log, gI, pts[1])
        @test size(gT) == size(ge)
        @test eltype(gT) == eltype(ge)

        for T in types
            gpts = convert.(T, pts)
            Xgpts = convert.(T, Xpts)
            @test compose(G, gpts[1], gpts[2]) ≈ gpts[1] * gpts[2]
            @test translate_diff(G, gpts[2], gpts[1], Xgpts[1], LeftAction()) ≈ Xgpts[1]
            @test translate_diff(G, gpts[2], gpts[1], Xgpts[1], RightAction()) ≈
                  transpose(gpts[2]) * Xgpts[1] * gpts[2]
            test_group(
                G,
                gpts,
                Xgpts,
                Xgpts;
                test_diff=true,
                test_invariance=true,
                test_lie_bracket=true,
                test_adjoint_action=true,
            )
        end

        @testset "Decorator forwards to group" begin
            DM = NotImplementedGroupDecorator(G)
            test_group(DM, pts, Xpts, Xpts; test_diff=true)
        end

        @testset "Group forwards to decorated" begin
            retraction_methods = [
                Manifolds.PolarRetraction(),
                Manifolds.QRRetraction(),
                Manifolds.GroupExponentialRetraction(LeftAction()),
                Manifolds.GroupExponentialRetraction(RightAction()),
            ]

            inverse_retraction_methods = [
                Manifolds.PolarInverseRetraction(),
                Manifolds.QRInverseRetraction(),
                Manifolds.GroupLogarithmicInverseRetraction(LeftAction()),
                Manifolds.GroupLogarithmicInverseRetraction(RightAction()),
            ]

            test_manifold(
                G,
                pts;
                test_injectivity_radius=false,
                test_project_tangent=true,
                test_musical_isomorphisms=false,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                exp_log_atol_multiplier=20,
                retraction_atol_multiplier=12,
                is_tangent_atol_multiplier=1.2,
                test_atlases=(Manifolds.RetractionAtlas(),),
            )

            @test injectivity_radius(G) == injectivity_radius(M)
            @test injectivity_radius(G, pts[1]) == injectivity_radius(M, pts[1])
            @test injectivity_radius(G, pts[1], PolarRetraction()) ==
                  injectivity_radius(M, pts[1], PolarRetraction())

            y = allocate(pts[1])
            exp!(G, y, pts[1], Xpts[1])
            @test isapprox(M, y, exp(M, pts[1], Xpts[1]))

            y = allocate(pts[1])
            retract!(G, y, pts[1], Xpts[1])
            @test isapprox(M, y, retract(M, pts[1], Xpts[1]))

            w = allocate(Xpts[1])
            inverse_retract!(G, w, pts[1], pts[2])
            @test isapprox(M, pts[1], w, inverse_retract(M, pts[1], pts[2]))

            w = allocate(Xpts[1])
            inverse_retract!(G, w, pts[1], pts[2], QRInverseRetraction())
            @test isapprox(
                M,
                pts[1],
                w,
                inverse_retract(M, pts[1], pts[2], QRInverseRetraction()),
            )

            if n == 2
                z = collect(reshape(1.0:4.0, 2, 2))
            else
                z = collect(reshape(1.0:9.0, 3, 3))
            end

            @test isapprox(M, project(G, z), project(M, z))
            z2 = allocate(pts[1])
            project!(G, z2, z)
            @test isapprox(M, z2, project(M, z))
        end

        @testset "vee/hat" begin
            X = Xpts[1]
            pe = Identity(G)

            Xⁱ = vee(G, identity_element(G), X)
            @test Xⁱ ≈ vee(G, pe, X)

            X2 = hat(G, pts[1], Xⁱ)
            @test isapprox(M, pe, X2, hat(G, pe, Xⁱ); atol=1e-6)
        end
        @testset "Identity and get_vector/get_coordinates" begin
            e = Identity(G)

            eF = Identity(SpecialEuclidean(n))
            if n == 2
                c = [1.0]
            else
                c = [1.0, 0.0, 0.0]
            end
            Y = zeros(representation_size(G))
            get_vector_lie!(G, Y, c, Manifolds.VeeOrthogonalBasis())
            @test Y ≈ get_vector_lie(G, c, Manifolds.VeeOrthogonalBasis())
            get_vector!(M, Y, identity_element(G), c, Manifolds.VeeOrthogonalBasis())
            @test Y ≈ get_vector(M, identity_element(G), c, Manifolds.VeeOrthogonalBasis())

            @test get_coordinates(
                M,
                identity_element(G),
                Y,
                Manifolds.VeeOrthogonalBasis(),
            ) == c
            c2 = similar(c)
            get_coordinates_lie!(G, c2, Y, Manifolds.VeeOrthogonalBasis())
            @test c == c2
            get_coordinates!(M, c2, identity_element(G), Y, Manifolds.VeeOrthogonalBasis())
            @test c == c2

            q = zeros(n, n)
            mul!(q, e, e)
            @test isone(q)
            e2 = Identity(G)
            @test mul!(e2, e, e) === e2
        end
    end
end


text:
module MadNLPSchur

import ..MadNLPGraph:
    @kwdef, Logger, @debug, @warn, @error,
    AbstractOptions, AbstractLinearSolver, EmptyLinearSolver, set_options!, SparseMatrixCSC, SubVector, StrideOneVector,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    introduce, factorize!, solve!, improve!, is_inertia, inertia,
    default_linear_solver, default_dense_solver, get_csc_view, get_cscsy_view, nnz, mul!,
    TwoStagePartition, set_blas_num_threads, blas_num_threads, @blas_safe_threads

const INPUT_MATRIX_TYPE = :csc

@kwdef mutable struct Options <: AbstractOptions
    schur_custom_partition::Bool = false
    schur_num_parts::Int = 2
    schur_part::Vector{Int} = Int[]
    schur_subproblem_solver::Module
    schur_dense_solver::Module
end

mutable struct SolverWorker
    V::Vector{Int}
    V_0_nz::Vector{Int}
    csc::SparseMatrixCSC{Float64,Int32}
    csc_view::SubVector{Float64}
    compl::SparseMatrixCSC{Float64,Int32}
    compl_view::SubVector{Float64}
    M::AbstractLinearSolver
    w::Vector{Float64}
end

mutable struct Solver <: AbstractLinearSolver
    csc::SparseMatrixCSC{Float64,Int32}
    inds::Vector{Int}
    tsp::TwoStagePartition

    schur::Matrix{Float64}
    colors
    fact

    V_0::Vector{Int}
    csc_0::SparseMatrixCSC{Float64,Int32}
    csc_0_view::SubVector{Float64}
    w_0::Vector{Float64}

    sws::Vector{SolverWorker}

    opt::Options
    logger::Logger
end

function Solver(csc::SparseMatrixCSC{Float64};
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                opt=Options(schur_subproblem_solver=default_linear_solver(),
                            schur_dense_solver=default_dense_solver()),
                logger=Logger(),
                kwargs...)

    set_options!(opt,option_dict,kwargs...)
    if string(opt.schur_subproblem_solver) == "MadNLP.Mumps"
        @warn(logger,"When Mumps is used as a subproblem solver, Schur is run in serial.")
        @warn(logger,"To use parallelized Schur, use Ma27 or Ma57.")
    end

    inds = collect(1:nnz(csc))
    tsp = TwoStagePartition(csc,opt.schur_part,opt.schur_num_parts)

    V_0   = findall(tsp.part.==0)
    colors = get_colors(length(V_0),opt.schur_num_parts)

    csc_0,csc_0_view = get_cscsy_view(csc,V_0,inds=inds)
    schur = Matrix{Float64}(undef,length(V_0),length(V_0))

    w_0 = Vector{Float64}(undef,length(V_0))

    sws = Vector{Any}(undef,opt.schur_num_parts)

    copied_option_dict = copy(option_dict)
    @blas_safe_threads for k=1:opt.schur_num_parts
        sws[k] = SolverWorker(
            tsp,V_0,csc,inds,k,opt.schur_subproblem_solver,logger,k==1 ? option_dict : copy(copied_option_dict))
    end
    fact = opt.schur_dense_solver.Solver(schur)
    return Solver(csc,inds,tsp,schur,colors,fact,V_0,csc_0,csc_0_view,w_0,sws,opt,logger)
end

get_colors(n0,K) = [findall((x)->mod(x-1,K)+1==k,1:n0) for k=1:K]

function SolverWorker(tsp,V_0,csc::SparseMatrixCSC{Float64},inds::Vector{Int},k,
                      SubproblemSolverModule::Module,logger::Logger,option_dict::Dict{Symbol,Any})

    V    = findall(tsp.part.==k)

    csc_k,csc_k_view = get_cscsy_view(csc,V,inds=inds)
    compl,compl_view = get_csc_view(csc,V,V_0,inds=inds)
    V_0_nz = findnz(compl.colptr)

    M    = length(V) == 0 ?
        EmptyLinearSolver() : SubproblemSolverModule.Solver(csc_k;option_dict=option_dict,logger=logger)
    w    = Vector{Float64}(undef,csc_k.n)

    return SolverWorker(V,V_0_nz,csc_k,csc_k_view,compl,compl_view,M,w)
end

function findnz(colptr)
    nz = Int[]
    for j=1:length(colptr)-1
        colptr[j]==colptr[j+1] || push!(nz,j)
    end
    return nz
end

function factorize!(M::Solver)
    M.schur.=0.
    M.csc_0.nzval.=M.csc_0_view
    M.schur.=M.csc_0
    @blas_safe_threads for sw in M.sws
        sw.csc.nzval.=sw.csc_view
        sw.compl.nzval.=sw.compl_view
        factorize!(sw.M)
    end

    # asymchronous multithreading doesn't work here
    for q = 1:length(M.colors)
        @blas_safe_threads for k = 1:length(M.sws)
            for j = M.colors[mod(q+k-1,length(M.sws))+1] # each subprob works on a different color
                factorize_worker!(j,M.sws[k],M.schur)
            end
        end
    end
    factorize!(M.fact)

    return M
end

function factorize_worker!(j,sw,schur)
    j in sw.V_0_nz || return
    sw.w.= view(sw.compl,:,j)
    solve!(sw.M,sw.w)
    mul!(view(schur,:,j),sw.compl',sw.w,-1.,1.)
end

function solve!(M::Solver,x::AbstractVector{Float64})
    M.w_0 .= view(x,M.V_0)
    @blas_safe_threads for sw in M.sws
        sw.w.=view(x,sw.V)
        solve!(sw.M,sw.w)
    end
    for sw in M.sws
        mul!(M.w_0,sw.compl',sw.w,-1.,1.)
    end
    solve!(M.fact,M.w_0)
    view(x,M.V_0).=M.w_0
    @blas_safe_threads for sw in M.sws
        x_view = view(x,sw.V)
        sw.w.= x_view
        mul!(sw.w,sw.compl,M.w_0,1.,1.)
        solve!(sw.M,sw.w)
        x_view.=sw.w
    end
    return x
end

is_inertia(M::Solver) = is_inertia(M.fact) && is_inertia(M.sws[1].M)
function inertia(M::Solver)
    numpos,numzero,numneg = inertia(M.fact)
    for k=1:M.opt.schur_num_parts
        _numpos,_numzero,_numneg =  inertia(M.sws[k].M)
        numpos += _numpos
        numzero += _numzero
        numneg += _numneg
    end
    return (numpos,numzero,numneg)
end

function improve!(M::Solver)
    for sw in M.sws
        improve!(sw.M) || return false
    end
    return true
end

function introduce(M::Solver)
    for sw in M.sws
        sw.M isa EmptyLinearSolver || return "schur equipped with "*introduce(sw.M)
    end
end

end # module


text:
include("bitstrings/mutations.jl")
include("bitstrings/initialpoint.jl")

include("graphs/utils.jl")

max_stars_repo_name: DTUComputeTONIA/SDRLS


text:
# ---
# title: 503. Next Greater Element II
# id: problem503
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/next-greater-element-ii/description/>
# hidden: true
# ---
#
# Given a circular array (the next element of the last element is the first
# element of the array), print the Next Greater Number for every element. The
# Next Greater Number of a number x is the first greater number to its
# traversing-order next in the array, which means you could search circularly to
# find its next greater number. If it doesn't exist, output -1 for this number.
#
# **Example 1:**
#
#
#
#     Input: [1,2,1]
#     Output: [2,-1,2]
#     Explanation: The first 1's next greater number is 2; The number 2 can't find next greater number; The second 1's next greater number needs to search circularly, which is also 2.
#
#
# **Note:** The length of given array won't exceed 10000.
#
#
## @lc code=start
using LeetCode

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


text:
Cassette.@context RecordingCtx
new_recording_ctx() = Cassette.disablehooks(RecordingCtx(metadata=AllocationRecord()))

function record_alloc!(record::AllocationRecord, val)
    push!(record.initial_sizes, size(val))
    push!(record.allocations, val)
end
record_alloc!(ctx::RecordingCtx, val) = record_alloc!(ctx.metadata, val)

@inline function Cassette.overdub(
    ctx::RecordingCtx, ::Type{Array{T,N}}, ::UndefInitializer, dims
) where {T,N}
    ret = Array{T, N}(undef, dims)
    record_alloc!(ctx, ret)
    return ret
end

function record_allocations(f, args...; kwargs...)
    ctx = new_recording_ctx()
    value = Cassette.overdub(ctx, f, args...; kwargs...)

    return (value=value, allocation_record=ctx.metadata,)
end


text:
module TestOneHot

using BangBang.Extras: broadcast_inplace!!
using InitialValues
using MicroCollections
using Test

@testset begin
    @test broadcast_inplace!!(+, InitialValue(+), OneHotArray(2 => 1, 3)) == [0, 1, 0]
    @test broadcast_inplace!!(+, [1, 2, 3], OneHotArray(2 => 1, 3)) == [1, 3, 3]
    @test broadcast_inplace!!(
        +,
        InitialValue(+),
        OneHotArray(2 => 1, 3),
        OneHotArray(3 => -2, 3),
    ) == [0, 1, -2]

    @test broadcast_inplace!!(*, InitialValue(*), OneHotArray(2 => 111, 3)) == [1, 111, 1]
    @test broadcast_inplace!!(*, [1, 2, 3], OneHotArray(2 => 111, 3)) == [1, 222, 3]
end

end  # module


text:
if normpath(@__DIR__) ∉ LOAD_PATH
    pushfirst!(LOAD_PATH, normpath(@__DIR__, "../.."))
    pushfirst!(LOAD_PATH, normpath(@__DIR__))
end

using DECAES
using DECAES.NormalHermiteSplines
using StaticArrays
using Random
using LinearAlgebra
const nhs = DECAES.NormalHermiteSplines

# Packages from this local env
using BenchmarkTools
using LaTeXStrings
using CairoMakie
set_theme!(theme_ggplot2(); resolution = (800,600), font = "CMU Serif")

function plot_neighbours(::Val{D} = Val(2)) where {D}
    grid  = DECAES.meshgrid(SVector{D,Float64}, [range(0, 1; length = 25) for _ in 1:D]...)
    surr  = DECAES.HermiteSplineSurrogate(I -> (1.0, zero(SVector{D,Float64})), grid)
    state = DECAES.DiscreteSurrogateSearcher(surr; mineval = 5, maxeval = typemax(Int))

    second(x) = x[2]
    xy(xs) = D == 1 ? (first.(xs), ones(length(xs))) : vec.((first.(xs), second.(xs)))

    while true
        xnew = rand(SVector{D,Float64})
        box = DECAES.minimal_bounding_box(state, xnew)
        oldseen = copy(state.seen)
        DECAES.evaluate_box!(surr, state, box; maxeval = typemax(Int))
        newseen = state.seen .!= oldseen

        fig = Figure()
        ax = fig[1,1] = Axis(fig)
        scatter!(xy(state.grid)...; markersize = 5, color = :black, axis = (; aspect = AxisAspect(1.0)), label = "grid")
        scatter!(xy(state.grid[oldseen])...; markersize = 20, color = :darkblue, label = "old pts")
        scatter!(xy(state.grid[newseen])...; markersize = 20, color = :orange, label = "new pts")
        scatter!(xy([xnew])...; markersize = 25, color = :red, marker = :star5, label = "query")
        fig[1,2] = Legend(fig, ax)
        display(fig)
        readline()
    end
end

function plot_bisection_search(
        ::Val{D},
        ::Val{ETL},
        prob::DECAES.NNLSDiscreteSurrogateSearch{D,T},
        opts::DECAES.T2mapOptions{T};
        surrtype = :hermite,
        npts     = opts.nRefAngles,
        mineval  = opts.nRefAnglesMin,
        maxeval  = npts,
    ) where {D,T,ETL}
    surrtype === :cubic && @assert D == 1 "cubic splines only support 1D"

    # build surrogate
    surr = surrtype === :cubic ?
        DECAES.CubicSplineSurrogate(prob) :
        DECAES.HermiteSplineSurrogate(prob)
    f_true = function(x...)
        α, β = length(x) == 1 ? (x[1], opts.SetRefConAngle) : (x[1], x[2])
        θ = DECAES.EPGOptions((; α = α, TE = opts.TE, T2 = 0.0, T1 = opts.T1, β = β), Val(ETL), T)
        A = DECAES.epg_decay_basis(θ, DECAES.logrange(opts.T2Range..., opts.nT2))
        nnls_prob = DECAES.NNLSProblem(A, prob.b)
        DECAES.solve!(nnls_prob, A, prob.b)
        return log(DECAES.chi2(nnls_prob))
    end

    # solve discrete search problem
    state = DECAES.DiscreteSurrogateSearcher(surr; mineval = mineval, maxeval = maxeval)
    minx, miny = DECAES.bisection_search(surr, state; maxeval = maxeval)
    # minx, miny = DECAES.local_search(surr, minx)

    # reconstruct surrogate from evaluated points and plot
    if surrtype === :cubic
        nodes = () -> surr.p
        values = () -> surr.u
        spl = DECAES.make_spline(first.(nodes()), values())
    else
        nodes = () -> map(x -> nhs._unnormalize(surr.spl, x), nhs._get_nodes(surr.spl))
        values = () -> nhs._get_values(surr.spl)
        spl = (x...) -> nhs.evaluate(surr.spl, SVector(x...))
    end

    if D == 1
        fig = Figure()
        ax = fig[1,1] = Axis(fig)
        lines!(surr.grid[1][1]..surr.grid[end][1], f_true; color = :darkblue, label = "f")
        lines!(surr.grid[1][1]..surr.grid[end][1], x -> spl(x); color = :darkred, label = "spl")
        scatter!(first.(nodes()), values(); markersize = 10, color = :blue, label = "pts")
        scatter!(first.(minx), [miny]; markersize = 10, color = :red, marker = :star5, label = "min")
        fig[1,2] = Legend(fig, ax)
        return fig
    else
        xs = LinRange(surr.grid[1,1][1], surr.grid[end,end][1], 100)
        ys = LinRange(surr.grid[1,1][2], surr.grid[end,end][2], 100)
        zs = spl.(xs, ys')
        zs_true = f_true.(xs, ys')
        zmin_true, Imin_true = findmin(zs_true)

        fig = Figure()
        ax1 = Axis(fig[1,1]; title = L"$\sigma(x)$")
        contourf!(ax1, xs, ys, zs; limits = extrema(zs_true), levels = 50)
        scatter!(ax1, (p->p[1]).(nodes()), (p->p[2]).(nodes()), values(); marker = :rect, markersize = 10, color = :black, label = "pts")
        scatter!(ax1, [minx[1]], [minx[2]]; markersize = 15, color = :red, marker = :star5, label = "spl min")

        ax2 = Axis(fig[1,2]; title = L"f(x)")
        cont2 = contourf!(ax2, xs, ys, zs_true; limits = extrema(zs_true), levels = 50)
        scatter!(ax2, (p->p[1]).(nodes()), (p->p[2]).(nodes()), values(); marker = :rect, markersize = 10, color = :black, label = "pts")
        scatter!(ax2, [xs[Imin_true[1]]], [ys[Imin_true[2]]]; markersize = 15, color = :white, marker = :star5, label = "true min")
        scatter!(ax2, [minx[1]], [minx[2]]; markersize = 15, color = :red, marker = :star5, label = "spl min")

        # Colorbar(fig[2,1:2], cont2; vertical = false, flipaxis = false)
        hidedecorations!.([ax1, ax2])
        return fig
    end
end
function plot_bisection_search(
        ::Val{D};
        flip   = 175.0,
        refcon = 150.0,
        npts   = 32,
        kwargs...,
    ) where {D}
    opts = DECAES.mock_t2map_opts(; MatrixSize = (1,1,1), nTE = 32, SetFlipAngle = flip, SetRefConAngle = refcon, nRefAngles = npts)
    prob = DECAES.mock_surrogate_search_problem(Val(D), Val(32), opts)
    plot_bisection_search(Val(D), Val(32), prob, opts; npts = npts, kwargs...)
end
function plot_bisection_search(b::AbstractVector, opts::T2mapOptions, ::Val{D}; kwargs...) where {D}
    prob = DECAES.mock_surrogate_search_problem(b, opts, Val(D))
    plot_bisection_search(Val(D), Val(opts.nTE), prob, opts; npts = opts.nRefAngles, kwargs...)
end

function benchmark_spline()
    x = sort(randn(MersenneTwister(0), 10))
    u = randn(MersenneTwister(0), 10)
    du = randn(MersenneTwister(0), 10)
    xi = Ref(mean(x))
    spl = nhs.interpolate(x, u, nhs.RK_H1())
    dspl = nhs.interpolate(x, u, x, du, nhs.RK_H1())
    @btime $(nhs._gram!)($(spl._gram.data), $(spl._nodes), $(spl._kernel))
    @btime $(nhs._gram!)($(dspl._gram.data), $(dspl._nodes), $(dspl._d_nodes), $(dspl._d_dirs), $(dspl._kernel))
    @btime $(nhs.evaluate)($spl, $xi[])
    @btime $(nhs.evaluate)($dspl, $xi[])
    @btime $(nhs.evaluate_derivative)($spl, $xi[])
    @btime $(nhs.evaluate_derivative)($dspl, $xi[])
end


text:
module Ex5_6
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.6 - Solution of an AFDP
println("Example 5.6")

# define s as an improper transfer function
s = rtf('s');
# define Gu(s), Gf(s) and Gw(s)
Gu = [(s+1)/(s+2); (s+2)/(s+3)];     # enter Gu(s)
Gw = [1/(s+2); 0];                   # enter Gw(s)
Gf = [(s+1)/(s+2) 0; 0 1];           # enter Gf(s)
p = 2; mu = 1; mw = 1; mf = 2;       # enter dimensions

Q = dss([eye(p) -Gu]); Rf = dss(Gf); Rw = dss(Gw);
atol = 1.e-7;

# check solvability using a random frequency
if minimum(maximum(abs.(evalfr(Rf,rand())),dims=1)) > 0.01

    # compress Rw to a full row rank matrix
    rw = rank(evalfr(Rw,rand())); nb = size(Q,1);
    if rw < nb
       h = ones(rw,nb);  # usable only for rw = 1
       # use alternatively h = rand(rw,nb);
       Q = h*Q; Rf = h*Rf; Rw = h*Rw;
    end

    # compute the quasi-co-outer-co-inner factorization
    Rwi, Rwo, info = goifac(Rw; atol)

    # compute optimal filter (for standard case)
    Q = gir(Rwo\Q; atol )              # update Q
    Rf = gir(Rwo\Rf; atol); Rw = Rwi;  # update Rf and Rw

    # check for poles on the extended imaginary axis
    poles = gpole([Q Rf])
    if any(isinf.(poles)) || minimum(abs.(real(poles))) < 0.0001
        # compute a stable and proper left coprime factorization
        # of [Q Rf Rw] with desired stability degree -3
        # update if solution is improper or unstable
        Q_Rf_Rw, M = glcf(gir([Q Rf Rw]; atol); atol, sdeg = -3, smarg = -3);

        # adjust denominator to unit infinity norm to match example
        Mnorm = ghinfnorm(M)[1]
        Q_Rf_Rw = Q_Rf_Rw/Mnorm
        Q = Q_Rf_Rw[:,1:p+mu]
        Rf = Q_Rf_Rw[:,p+mu+1:p+mu+mf]
        Rw = Q_Rf_Rw[:,p+mu+mf+1:end]
        println("Q = $(dss2rm(Q, atol = 1.e-7))")
        println("Rf = $(dss2rm(Rf, atol = 1.e-7))")
        println("Rw = $(dss2rm(Rw, atol = 1.e-7))")
        gap = fdhinfminus(Rf)[1]
        println("gap = $gap")
    end
else
    @info "No solution exists"
end

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
"""
 Convert cartesian coordinates to spherical coordinates (θ,ϕ,r)
 θ, ϕ are in radians
 ϕ is the latitude, i.e., the angle from the equator.
"""
function xyz2θϕr(x)
  r = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
  θ = atan(x[2], x[1])
  ϕ = asin(x[3]/r)
  VectorValue(θ,ϕ,r)
end

"""
  Map point on the unit sphere on ambientspace to the parametric space of the unit sphere
  (x,y,z) ↦ (θ,ϕ)
"""
function xyz2θϕ(x)
  θ = atan(x[2], x[1])
  ϕ = asin(x[3])
  VectorValue(θ,ϕ)
end

"""
  Map point on the parametric space of the unit sphere to the ambient space
  (θ,ϕ) ↦ (x,y,z)
"""
function θϕ2xyz(θϕ)
  θ,ϕ = θϕ
  x = cos(θ)*cos(ϕ)
  y = sin(θ)*cos(ϕ)
  z = sin(ϕ)
  VectorValue(x,y,z)
end

"""
Matrix transformation from spherical vector field to Cartesian vector field.
  θ∈(0,2π)
  ϕ∈(-π/2,π/2)
  r=constant for the sphere
"""
function spherical_to_cartesian_matrix(θϕr)
  θ,ϕ,r = θϕr
  TensorValue(-sin(θ)       , cos(θ)       ,      0,
              -sin(ϕ)*cos(θ),-sin(ϕ)*sin(θ), cos(ϕ),
               cos(ϕ)*cos(θ), cos(ϕ)*sin(θ), sin(ϕ))
end

"""
Matrix transformation from cartesian vector field to spherical vector field
"""
function cartesian_to_spherical_matrix(xyz)
  x,y,z = xyz
  sr = sqrt(x^2+y^2+z^2)
  cr = sqrt(x^2+y^2)
  TensorValue(-y/cr, (x*z)/(sr*cr), x/sr,
               x/cr, (y*z)/(sr*cr), y/sr,
                  0,   -cr/(sr*cr), z/sr)
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function GitIndex(repo::GitRepo)
    idx_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_repository_index, :libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), idx_ptr_ptr, repo.ptr)
    return GitIndex(idx_ptr_ptr[])
end

function read!(idx::GitIndex, force::Bool = false)
    @check ccall((:git_index_read, :libgit2), Cint, (Ptr{Void}, Cint), idx.ptr, Cint(force))
    return idx
end

function write!(idx::GitIndex)
    @check ccall((:git_index_write, :libgit2), Cint, (Ptr{Void},), idx.ptr)
    return idx
end

function write_tree!(idx::GitIndex)
    oid_ptr = Ref(Oid())
    @check ccall((:git_index_write_tree, :libgit2), Cint,
                 (Ptr{Oid}, Ptr{Void}), oid_ptr, idx.ptr)
    return oid_ptr[]
end

function owner(idx::GitIndex)
    repo_ptr = ccall((:git_index_owner, :libgit2), Ptr{Void},
                     (Ptr{Void},), idx.ptr)
    return GitRepo(repo_ptr)
end

function read_tree!(idx::GitIndex, tree_id::Oid)
    repo = owner(idx)
    tree = get(GitTree, repo, tree_id)
    try
        @check ccall((:git_index_read_tree, :libgit2), Cint,
                     (Ptr{Void}, Ptr{Void}), idx.ptr, tree.ptr)
    finally
        finalize(tree)
    end
end

function add!{T<:AbstractString}(idx::GitIndex, files::T...;
             flags::Cuint = Consts.INDEX_ADD_DEFAULT)
    sa = StrArrayStruct(files...)
    try
        @check ccall((:git_index_add_all, :libgit2), Cint,
                     (Ptr{Void}, Ptr{StrArrayStruct}, Cuint, Ptr{Void}, Ptr{Void}),
                      idx.ptr, Ref(sa), flags, C_NULL, C_NULL)
    finally
        finalize(sa)
    end
end

function update!{T<:AbstractString}(idx::GitIndex, files::T...)
    sa = StrArrayStruct(files...)
    try
        @check ccall((:git_index_update_all, :libgit2), Cint,
                     (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{Void}, Ptr{Void}),
                      idx.ptr, Ref(sa), C_NULL, C_NULL)
    finally
        finalize(sa)
    end
end

function remove!{T<:AbstractString}(idx::GitIndex, files::T...)
    sa = StrArrayStruct(files...)
    try
        @check ccall((:git_index_remove_all, :libgit2), Cint,
                     (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{Void}, Ptr{Void}),
                      idx.ptr, Ref(sa), C_NULL, C_NULL)
    finally
        finalize(sa)
    end
end

function add!{T<:AbstractString}(repo::GitRepo, files::T...;
             flags::Cuint = Consts.INDEX_ADD_DEFAULT)
    with(GitIndex, repo) do idx
        add!(idx, files..., flags = flags)
        write!(idx)
    end
    return
end

function update!{T<:AbstractString}(repo::GitRepo, files::T...)
    with(GitIndex, repo) do idx
        update!(idx, files...)
        write!(idx)
    end
    return
end

function remove!{T<:AbstractString}(repo::GitRepo, files::T...)
    with(GitIndex, repo) do idx
        remove!(idx, files...)
        write!(idx)
    end
    return
end

function read!(repo::GitRepo, force::Bool = false)
    with(GitIndex, repo) do idx
        read!(idx, force)
    end
    return
end

function Base.count(idx::GitIndex)
    return ccall((:git_index_entrycount, :libgit2), Csize_t, (Ptr{Void},), idx.ptr)
end

function Base.getindex(idx::GitIndex, i::Csize_t)
    ie_ptr = ccall((:git_index_get_byindex, :libgit2), Ptr{Void},
                  (Ptr{Void}, Csize_t), idx.ptr, i-1)
    ie_ptr == C_NULL && return nothing
    return unsafe_load(convert(Ptr{IndexEntry}, ie_ptr), 1)
end
Base.getindex(idx::GitIndex, i::Int) = getindex(idx, Csize_t(i))

function Base.find(path::String, idx::GitIndex)
    pos_ref = Ref{Csize_t}(0)
    ret = ccall((:git_index_find, :libgit2), Cint,
                  (Ref{Csize_t}, Ptr{Void}, Cstring), pos_ref, idx.ptr, path)
    ret == Error.ENOTFOUND && return Nullable{Csize_t}()
    return Nullable(pos_ref[]+1)
end

stage(ie::IndexEntry) = ccall((:git_index_entry_stage, :libgit2), Cint, (Ptr{IndexEntry},), Ref(ie))

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
export Arc

mobius(a,b,c,d,z) = (a*z+b)./(c*z+d)
mobius(a,b,c,d,z::Number) = isinf(z)?a/c:(a*z+b)./(c*z+d)
mobiusinv(a,b,c,d,z) = mobius(d,-b,-c,a,z)

mobiusD(a,b,c,d,z) = (a*d-b*c)./(d+c*z).^2
mobiusinvD(a,b,c,d,z) = mobiusD(d,-b,-c,a,z)

immutable Arc{T,V<:Real,TT} <: IntervalDomain{TT}
    center::T
    radius::V
    angles::Tuple{V,V}
    Arc(c,r,a) = new(T(c),V(r),Tuple{V,V}(a))
end

Arc{T<:Number,V<:Real,V1<:Real,V2<:Real}(c::T,r::V,t::Tuple{V1,V2}) =
    Arc{promote_type(T,V,V1,V2),
        promote_type(real(T),V,V1,V2),
        Complex{promote_type(real(T),V,V1,V2)}}(c,r,t)
Arc{T<:Number,V<:Real}(c::Vec{2,T},r::V,t::Tuple{V,V}) =
    Arc{Vec{2,promote_type(T,V)},
        promote_type(real(T),V),
        Vec{2,promote_type(real(T),V)}}(c,r,t)
Arc{T<:Number,V<:Real,V1<:Real,V2<:Real}(c::Vec{2,T},r::V,t::Tuple{V1,V2}) =
    Arc{Vec{2,promote_type(T,V,V1,V2)},
        promote_type(real(T),V,V1,V2),
        Vec{2,promote_type(T,V,V1,V2)}}(c,r,t)

Arc(c::Tuple,r,t) = Arc(Vec(c...),r,t)
Arc(c,r,t0,t1) = Arc(c,r,(t0,t1))

Base.complex{V<:Vec}(a::Arc{V}) = Arc(complex(a.center...),a.radius,a.angles)

isambiguous(d::Arc) =
    isnan(d.center) && isnan(d.radius) && isnan(d.angles[1]) && isnan(d.angles[2])
Base.convert{T<:Number,V<:Real}(::Type{Arc{T,V}},::AnyDomain) =
    Arc{T,V}(NaN,NaN,(NaN,NaN))
Base.convert{IT<:Arc}(::Type{IT},::AnyDomain) =
    Arc(NaN,NaN,(NaN,NaN))

arclength(d::Arc) = d.radius*(d.angles[2]-d.angles[1])

function mobiuspars(z0,r,t0,t1)
    c=exp(im*t0/2)+exp(im*t1/2)
    f=exp(im*t0/2)-exp(im*t1/2)
    d=exp(im*t0/2+im*t1/2)*r
    -c,c*(z0+d),f,f*(d-z0)
end

mobiuspars(a::Arc) = mobiuspars(a.center,a.radius,a.angles...)

for OP in (:mobius,:mobiusinv,:mobiusD,:mobiusinvD)
    @eval $OP(a::Arc,z) = $OP(mobiuspars(a)...,z)
end

tocanonical{T<:Number}(a::Arc{T},x) = real(mobius(a,x))
tocanonicalD{T<:Number}(a::Arc{T},x) = mobiusD(a,x)
fromcanonical{T<:Number,V<:Real,TT<:Complex}(d::Arc{T,V,TT},v::AbstractArray) =
    eltype(d)[fromcanonical(d,vk) for vk in v]
fromcanonical{T<:Number,V<:Real,TT<:Complex}(a::Arc{T,V,TT},x) =
    mobiusinv(a,x)
fromcanonicalD{T<:Number}(a::Arc{T},x) = mobiusinvD(a,x)

tocanonical{V<:Vec}(a::Arc{V},x::Vec) =
    tocanonical(complex(a),complex(x...))
fromcanonical{V<:Vec}(a::Arc{V},x::Number) =
    Vec(reim(fromcanonical(complex(a),x))...)
fromcanonicalD{V<:Vec}(a::Arc{V},x::Number) =
    Vec(reim(fromcanonicalD(complex(a),x))...)

## information

Base.issubset(a::Arc,b::Arc) =
    a.center==b.center && a.radius==b.radius &&
                    (b.angles[1]≤a.angles[1]≤a.angles[2]≤b.angles[2] ||
                            b.angles[1]≥a.angles[1]≥a.angles[2]≥b.angles[2])

# Algebra

for op in (:*,:.*)
    @eval begin
        $op(c::Real,d::Arc) =
            Arc(c*d.center,c*d.radius,(sign(c)*d.angles[1],sign(c)*d.angles[2]))
        $op(d::Arc,c::Real) = $op(c,d)
    end
end

for op in (:+,:-,:.+,:.-)
    @eval begin
        $op(c::Number,d::Arc) = Arc($op(c,d.center),d.radius,d.angles)
        $op(d::Arc,c::Number) = Arc($op(d.center,c),d.radius,d.angles)
    end
end

# allow exp(im*Interval(0,1)) for constructing arc
function Base.exp{CMP<:Complex}(d::Interval{CMP})
    @assert isapprox(real(d.a),0) && isapprox(real(d.b),0)
    Arc(0,1,(imag(d.a),imag(d.b)))
end


text:
@testset "JJ terms" begin
    @testset "_terms_jw" begin
        function terms_reference(j::HalfInteger, w::Integer)
            w <= 2j+1 || throw(ArgumentError("w=$w too large for $orb orbital"))

            2w ≥ 2j+1 && (w = convert(Int, 2j) + 1 - w)
            w == 0 && return [zero(HalfInt)]
            w == 1 && return [j]

            # Forms full Cartesian product of all mⱼ, not necessarily the most
            # performant method.
            mⱼs = filter(allunique, collect(AtomicLevels.allchoices([-j:j for i = 1:w])))
            MJs = map(x -> reduce(+, x), mⱼs) # TODO: make map(sum, mⱼs) work

            Js = HalfInt[]

            while !isempty(MJs)
                # Identify the maximum MJ and associate it with J.
                MJmax = maximum(MJs)
                N = count(isequal(MJmax), MJs)
                append!(Js, repeat([MJmax], N))
                # Remove the -MJ:MJ series, N times.
                for MJ = -MJmax:MJmax
                    deleteat!(MJs, findall(isequal(MJ), MJs)[1:N])
                end
            end

            # Do we really want unique here?
            sort(unique(Js))
        end

        for twoj = 0:10, w = 1:twoj+1
            j = HalfInteger(twoj//2)
            ts = AtomicLevels._terms_jw(j, w)

            @test issorted(ts, rev=true) # make sure that the array is sorted in descending order
            @test terms_reference(j, w) == sort(unique(ts))
            # Check particle-hole symmetry
            if w != twoj + 1
                @test AtomicLevels._terms_jw(j, twoj+1-w) == ts
            end
        end
    end

    @testset "jj coupling of equivalent electrons" begin
        @test terms(ro"1s", 0) == [0]
        @test terms(ro"1s", 1) == [1//2]
        @test terms(ro"1s", 2) == [0]

        @test terms(ro"3d-", 0) == [0]
        @test terms(ro"3d-", 1) == [3//2]
        @test terms(ro"3d-", 4) == [0]

        @test terms(ro"Xd", 0) == [0]
        @test terms(ro"Xd", 1) == [5//2]
        @test terms(ro"Xd", 6) == [0]

        # Table 4.5, Cowan 1981
        foreach([
            (ro"1s",ro"2p-") => [(0,2) => 0, 1 => 1//2],
            (ro"2p",ro"3d-") => [(0,4) => 0, (1,3) => 3//2, 2 => [0,2]],
            (ro"3d",ro"4f-") => [(0,6) => 0, (1,5) => 5//2, (2,4) => [0,2,4],
                                 3 => [3//2,5//2,9//2]],
            (ro"4f",ro"5g-") => [(0,8) => 0, (1,7) => 7//2, (2,6) => [0,2,4,6],
                                 (3,5) => [3//2,5//2,7//2,9//2,11/2,15//2],
                                 4 => [0,2,2,4,4,5,6,8]],
            (ro"5g",ro"6h-") => [(0,10) => 0, (1,9) => 9//2, (2,8) => [0,2,4,6,8],
                                 (3,7) => [3//2,5//2,7//2,9//2,9//2,11//2,13//2,15//2,17//2,21//2],
                                 (4,6) => [0,0,2,2,3,4,4,4,5,6,6,6,7,8,8,9,10,12],
                                 5 => [1//2,3//2,5//2,5//2,7//2,7//2,9//2,9//2,9//2,11//2,11//2,13//2,13//2,15//2,15//2,17//2,17//2,19//2,21//2,25//2]]
        ]) do (orbs,wsj)
            foreach(orbs) do orb
                foreach(wsj) do (ws,j)
                    js = map(HalfInteger, isa(j, Number) ? [j] : j)
                    foreach(ws) do w
                        @test terms(orb,w) == js
                    end
                end
            end
        end
    end
end


text:
# This file was generated, do not modify it. # hide
HotTreeClf = @pipeline(OneHotEncoder(),
                       DecisionTreeClassifier())

mdl = HotTreeClf
mach = machine(mdl, X, y)
fit!(mach);


text:
module GraphicsManager

using FileManager
using FileIO.query

export IGraphicsReference
export IGraphicsData
export IGraphicsShaderProgram
export IGraphicsShader
export IGraphicsShaderProperty

""" TODO """
abstract type IGraphicsReference end

""" TODO """
abstract type IGraphicsData <: IGraphicsReference end

""" TODO """
abstract type IGraphicsShaderProgram <: IGraphicsReference end

""" TODO """
abstract type IGraphicsShader <: IGraphicsReference end

""" TODO """
abstract type IGraphicsShaderProperty <: IGraphicsReference end

export AbstractGraphicsReference
export AbstractGraphicsData
export AbstractGraphicsShaderProgram
export AbstractGraphicsShader
export AbstractGraphicsShaderProperty

""" TODO """
AbstractGraphicsReference = Union{Void,IGraphicsReference}

""" TODO """
AbstractGraphicsData = Union{Void,IGraphicsData}

""" TODO """
AbstractGraphicsShaderProgram = Union{Void,IGraphicsShaderProgram}

""" TODO """
AbstractGraphicsShader = Union{Void,IGraphicsShader}

""" TODO """
AbstractGraphicsShaderProperty = Union{Void,IGraphicsShaderProperty}

""" TODO """
type Manager
	handle ::Any
	Manager() = new()
end

MANAGER = Manager()

""" TODO """
setGraphicsDriver(o::Any) = (MANAGER.handle = o)

""" TODO """
getGraphicsDriver() = MANAGER.handle

""" TODO """
GRAPHICSDRIVER() = getGraphicsDriver()

export GRAPHICSDRIVER

end #GraphicsManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
export readlistmat, readlistmat!, readevalues, writelistmat

"""
    readlistmat(fd::IO, nodes1::Vector, nodes2::Vector; <keyword arguments>)
    readlistmat(file::AbstractString, nodes1::Vector, nodes2::Vector; <keyword arguments>)

Reads a numerical matrix stored in list format, where the first and second columns correspond
to string vectors nodes1 and nodes2, respectively. E.g.

```
nodeA1 nodeA2 4.5
nodeB1 nodeB2 3.4
nodeA1 nodeB2 0.3
nodeB1 nodeA2 0.6
```

Returns a dense matrix by default. Set keyword option `dense=false` to return a sparse matrix.

# Arguments
- `fd`,`file` : file name or file I/O
- `nodes1`,`nodes1` : node vectors corresponding to 1st and 2nd columns

# Keyword arguments
- `header=false` : set to true to ignore first line
- `ignore=false` : set to true to ignore nodes in file that is not in nodes1 or nodes2
- `dense=true` : set to false to return sparse matrix
"""
function readlistmat(fd::IO, nodes1::Vector{<:AbstractString},
                     nodes2::Vector{<:AbstractString};
                     header=false, ignore=false, minval=-Inf,
                     dense=true,missing=0.0)
    if dense
        return readlistmat!(fd, fill(missing,length(nodes1),length(nodes2)),
                            nodes1, nodes2,
                            header=header, ignore=ignore,
                            minval=minval)
    end
    nodes1 = indexmap(nodes1)
    nodes2 = indexmap(nodes2)
    I = Int[]
    J = Int[]
    V = Float64[]
    header && readline(fd)
    for line in eachline(fd)
        vals = split(strip(line))
        if !ignore ||
            (haskey(nodes1,vals[1]) && haskey(nodes2,vals[2]))
            push!(I, nodes1[vals[1]])
            push!(J, nodes2[vals[2]])
            push!(V, max(minval, parse(Float64,vals[3])))
        end
    end
    sparse(I,J,V,length(nodes1),length(nodes2))
end
readlistmat(file::AbstractString,args...;vargs...) =
    open(fd -> readlistmat(fd,args...;vargs...), file, "r")

"""
    readlistmat!(fd::IO, B::AbstractMatrix, nodes1::Vector, nodes2::Vector; <keyword arguments>)
    readlistmat!(file::AbstractString, B::AbstractMatrix, nodes1::Vector, nodes2::Vector; <keyword arguments>)

Same as [`readlistmat`](@ref) but stores the result in B.
"""
function readlistmat!(fd::IO, B::AbstractMatrix, nodes1::Vector{<:AbstractString},
                      nodes2::Vector{<:AbstractString};
                      header=false, ignore=false, minval=-Inf)
    nodes1 = indexmap(nodes1)
    nodes2 = indexmap(nodes2)
    header && readline(fd)
    for line in eachline(fd)
        vals = split(strip(line))
        if !ignore || (haskey(nodes1,vals[1]) && haskey(nodes2,vals[2]))
            B[nodes1[vals[1]],nodes2[vals[2]]] = max(minval,parse(Float64,vals[3]))
        end
    end
    B
end
readlistmat!(file::AbstractString,args...;vargs...) =
    open(fd -> readlistmat!(fd,args...;vargs...), file, "r")

readevalues(fd::IO, args...; vargs...) = readlistmat(fd, args...; missing=Inf, vargs...)
readevalues(file::AbstractString,args...;vargs...) =
    open(fd -> readevalues(fd,args...;vargs...), file, "r")

function writelistmat(fd::IO, S::SparseMatrixCSC, nodes1::Vector{<:AbstractString},
                      nodes2::Vector{<:AbstractString})
    I,J,V = findnz(S)
    for (i,j,v) in zip(I,J,V)
        println(fd, "$(nodes1[i]) $(nodes2[j]) $v")
    end
end
writelistmat(file::AbstractString, args...;vargs...) =
    open(fd -> writelistmat(fd, args...;vargs...), file, "w")


text:
using Test
using Plots, Parameters, OrdinaryDiffEq

# hip extensor parameters, single ramp from 0.5 to 1.0 over 0.1 s
cc_p = CCParameters(838.2, 645.0, 26.0, 7.94, 1.0, 1.0, 1.0, 4.93, 1.64)
sec_p = SECParameters(2854.48)
α_p = ActivationProfile(0.5, ActivationRamp(0.0, 0.1, 1.0))

## TorqueGenerator
include("torque_generator.jl")
θ₀ = π / 2    # inital torque generator angle
ω₀ = 2π     # intial torque generator angular velocity
tqgen = TorqueGenerator(θ₀, ω₀, α_p, cc_p, sec_p)
@testset "torque generator" verbose = true begin
    # @test abs(tqgen.sec.τ - tqgen.cc.τ * tqgen.α) < 1e-6 # initial CC and SEC matching
    error = 0
    for θ₀ = 0:0.001:2π
        tqgen = TorqueGenerator(θ₀, ω₀, α_p, cc_p, sec_p)
        match = !(abs(tqgen.sec.τ - tqgen.cc.τ * tqgen.α) < 1e-6)
        error += match
    end
    @test error == 0
    θ₀ = rand(0.0:0.001:2π)
    tqgen = TorqueGenerator(θ₀, ω₀, α_p, cc_p, sec_p)
    @test tqgen.θ == tqgen.sec.θ + tqgen.cc.θ
    @test tqgen.ω == tqgen.sec.ω + tqgen.cc.ω
end

# test integration - he inital conditions
θ₀ = 215.89 |> deg2rad
ω₀ = -429.803 |> deg2rad
cc_p = CCParameters(838.2, 645.0, 26.0, 7.94, 1.0, 1.0, 1.0, 4.93, 1.64)
sec_p = SECParameters(2854.48)
α_p = ActivationProfile(0.9994277764, 0.6317965682e-01, 0.2115428829, 0.9106166714, 0.8071698914e-03, 0.1604915461, 0.1002353143)
tqgen = TorqueGenerator(θ₀, ω₀, cc_p, sec_p, α_p)
# θcc = tqgen.cc.θ
# t = 0.000
# dt = 0.001
# while t < 0.01
#     println("original values: t=$t, θ=$θ₀, ω=$ω₀, θcc=$θcc, ωcc=$(tqgen.cc.ω) τ=$(-tqgen.τ)")
#     τ = -update_torque_generator!(tqgen, t, θ₀, ω₀, θcc)
#     θ₀ += ω₀ * dt
#     ω₀ += τ * dt
#     θcc = tqgen.cc.θ + tqgen.cc.ω * dt
#     t += dt
#     println("updated values: t=$t, θ=$θ₀, ω=$ω₀, θcc=$θcc, ωcc=$(tqgen.cc.ω), τ=$(-tqgen.τ)")

# end

u₀ = [θ₀, ω₀, tqgen.cc.θ]
f(u, p, t) = begin
    τ = update_torque_generator!(p, t, u...)
    return [u[2], τ, p.cc.ω]
end
prob = ODEProblem(f, u₀, (0.0, 0.111), tqgen)
sol = solve(prob, Tsit5(), saveat = 0.001)
tq = [torque(sol.prob.p, sol(t)[1], sol(t)[3]) for t in sol.t]

## CC
include("cc.jl")
# plot relationships
plot(cc_p)

## activations
t1 = collect(0:0.001:0.2);
t2 = collect(0:0.001:0.1);
t3 = collect(0:0.001:0.2);
r1 = ramp.(t1, 0.0, 0.2, 1.0);
r2 = ramp.(t2, 1.0, 0.1, 0.5);
r3 = ramp.(t3, 0.5, 0.2, 0.0);
t_true = collect(0:0.001:1.0)
a_true = [repeat([0.0], 99); r1; repeat([1.0], 199); r2; r3; repeat([0.0], 199)]

a0 = 0.0
r1 = ActivationRamp(0.1, 0.2, 1.0)
r2 = ActivationRamp(0.2, 0.1, 0.5)
r3 = ActivationRamp(0.0, 0.2, 0.0)
p = ActivationProfile(a0, [r1, r2, r3])

## PEC
include("pec.jl")
# recreate figures from Reiner and Edrich (1999)
ϕh = range(-40, 120, length = 1000)
ϕk = range(-10, 150, length = 1000)
ϕa = range(-30, 60, length = 1000)

# hip
knee_angles = [0, 45, 90]
hip_plt = plot(title = "hip", xticks = [-40:20:120...], yticks = [-80:20:60...], ylims = (-90, 70))
for knee_angle in knee_angles
    plot!(ϕh, torqueHip.(ϕh, knee_angle), label = "ϕk = $knee_angle")

end
# knee
hip_angles = [0, 45, 90, 120]
knee_plt = plot(title = "knee", ylims = (-30, 30), xticks = [0:50:150...], yticks = [-20:10:30...])
for hip_angle in hip_angles
    plot!(ϕk, torqueKnee.(hip_angle, ϕk, 0), label = "ϕh = $hip_angle")
end
# ankle
knee_angles = [0, 60]
ankle_plt = plot(title = "ankle", ylims = (-20, 30), xticks = [-30:10:60...], yticks = [-20:5:30...])
for knee_angle in knee_angles
    plot!(ϕa, torqueAnkle.(knee_angle, ϕa), label = "ϕk = $knee_angle")
end
plot(hip_plt, knee_plt, ankle_plt)

θh = deg2rad.(0:0.1:220) # range(50, 220, length=1000) .|> deg2rad
θk = deg2rad.(0:0.1:220) # range(0, 180, length=1000) .|> deg2rad
θa = deg2rad.(0:0.1:220) # range(50, 180, length=1000) .|> deg2rad
# hip
knee_angles = [45, 90, 180] .|> deg2rad
hip_plt = plot(title = "hip", ylims = (-100, 100), legend = :topleft)
for knee_angle in knee_angles
    plot!(rad2deg.(θh), torqueHip.(convertHip.(θh), convertKnee(knee_angle)), label = "θk = $(rad2deg(knee_angle))")
end
hip_plt
# knee
hip_angles = [45, 90, 135, 180] .|> deg2rad
knee_plt = plot(title = "knee", ylims = (-100, 100), legend = :topleft)
for hip_angle in hip_angles
    plot!(rad2deg.(θk), torqueKnee.(convertHip(hip_angle), convertKnee.(θk), convertKnee.(0)), label = "θh = $(rad2deg(hip_angle))")
end
# ankle
knee_angles = [45, 135] .|> deg2rad
ankle_plt = plot(title = "ankle", ylims = (-100, 100))
for knee_angle in knee_angles
    plot!(rad2deg.(θa), torqueAnkle.(convertKnee(knee_angle), convertAnkle.(θa)), label = "θk = $(rad2deg(knee_angle))")
end
plot(hip_plt, knee_plt, ankle_plt)


text:
"""
    galerkin_tensor(
        B::AbstractBSplineBasis,
        (D₁::Derivative, D₂::Derivative, D₃::Derivative),
        [T = Float64],
    )

Compute 3D banded tensor appearing from quadratic terms in Galerkin method.

As with [`galerkin_matrix`](@ref), it is also possible to combine different
functional bases by passing, instead of `B`, a tuple `(B₁, B₂, B₃)` of three
`AbstractBSplineBasis`.
For now, the first two bases, `B₁` and `B₂`, must have the same length.

The tensor is efficiently stored in a [`BandedTensor3D`](@ref) object.
"""
function galerkin_tensor end

function galerkin_tensor(
        Bs::NTuple{3,AbstractBSplineBasis},
        deriv::DerivativeCombination{3},
        ::Type{T} = Float64,
    ) where {T}
    _check_bases(Bs)
    dims = length.(Bs)
    b = order(first(Bs)) - 1   # band width
    if length(Bs[1]) != length(Bs[2])
        throw(DimensionMismatch("the first two bases must have the same lengths"))
    end
    δl, δr = num_constraints(Bs[3]) .- num_constraints(Bs[1])
    A = BandedTensor3D{T}(undef, dims, Val(b), bandshift=(0, 0, δl))
    galerkin_tensor!(A, Bs, deriv)
end

galerkin_tensor(B::AbstractBSplineBasis, args...) =
    galerkin_tensor((B, B, B), args...)

"""
    galerkin_tensor!(
        A::BandedTensor3D,
        B::AbstractBSplineBasis,
        (D₁::Derivative, D₂::Derivative, D₃::Derivative),
    )

Compute 3D Galerkin tensor in-place.

See [`galerkin_tensor`](@ref) for details.
"""
function galerkin_tensor! end

function galerkin_tensor!(
        A::BandedTensor3D, Bs::NTuple{3,AbstractBSplineBasis},
        deriv::DerivativeCombination{3},
    )
    _check_bases(Bs)

    Ns = size(A)
    if Ns != length.(Bs)
        throw(DimensionMismatch("wrong dimensions of Galerkin tensor"))
    end

    Bi, Bj, Bl = Bs
    Ni, Nj, Nl = Ns
    @assert Ni == Nj  # verified at construction of BandedTensor3D

    # Orders and knots are assumed to be the same (see _check_bases).
    k = order(Bi)
    ts = knots(Bi)

    if bandwidth(A) != k - 1
        throw(ArgumentError("BandedTensor3D must have bandwidth = $(k - 1)"))
    end

    δl, _ = num_constraints(Bl) .- num_constraints(Bi)

    if bandshift(A) != (0, 0, δl)
        throw(ArgumentError("BandedTensor3D must have bandshift = (0, 0, $δl)"))
    end

    same_12 = Bi == Bj && deriv[1] == deriv[2]
    same_13 = Bi == Bl && deriv[1] == deriv[3]
    same_23 = Bj == Bl && deriv[2] == deriv[3]

    # Quadrature information (weights, nodes).
    quadx, quadw = _quadrature_prod(Val(3k - 3))

    fill!(A, 0)
    T = eltype(A)
    Al = @MMatrix zeros(T, 2k - 1, 2k - 1)
    nlast = last(eachindex(ts))

    @inbounds for n in eachindex(ts)
        n == nlast && break
        tn, tn1 = ts[n], ts[n + 1]
        tn1 == tn && continue  # interval of length = 0

        metric = QuadratureMetric(tn, tn1)
        xs = metric .* quadx

        is = nonzero_in_segment(Bi, n)
        js = nonzero_in_segment(Bj, n)
        ls = nonzero_in_segment(Bl, n)

        bis = eval_basis_functions(Bi, is, xs, deriv[1])
        bjs = same_12 ?
            bis : eval_basis_functions(Bj, js, xs, deriv[2])
        bls = same_13 ?
            bis : same_23 ?
            bjs : eval_basis_functions(Bl, ls, xs, deriv[3])

        # We compute the submatrix A[:, :, l] for each `l`.
        for (nl, l) in enumerate(ls)
            fill!(Al, 0)
            inds = BandedTensors.band_indices(A, l)
            # @assert is ⊆ inds && js ⊆ inds
            # @assert size(Al, 1) == size(Al, 2) == length(inds)
            δi = searchsortedlast(inds, is[1]) - 1
            δj = searchsortedlast(inds, js[1]) - 1
            for nj in eachindex(js), ni in eachindex(is)
                Al[ni + δi, nj + δj] =
                    metric.α * ((bis[ni] .* bjs[nj] .* bls[nl]) ⋅ quadw)
            end
            A[:, :, l] += Al
        end
    end

    A
end


text:
using JSON

#Constants
const e_charge = 1.602176453E-09
const e_mass = 9.1093826E-31
const c_0 = 299792458
const ε_0 = 8.85419E-12
const μ_0 = 4E-07

#Reading the dielectrics files that defines de problem
dielectrics = JSON.Parser.parsefile(".\\src\\Dielectrics.json")

#Reading the config files that defines de problem
scheme_prob = JSON.Parser.parsefile(".\\src\\Scheme.json")

for i in range(1, 20, step=1)
    println(i)
end