text:
module CommandLine

# TODO: accept command-line arguments to set the instance parsing configuration

# External packages used.
using UnPack: @unpack
import TimerOutputs.@timeit
using JuMP
import MathOptInterface
const MOI = MathOptInterface
import Dates
import Dates.@dateformat_str
import ArgParse
import ArgParse: ArgParseSettings, set_default_arg_group!, add_arg_group!
import ArgParse: @add_arg_table!, add_arg_table!

include("./SolversArgs.jl") # empty_configured_model, *parse_settings()
using .SolversArgs

import ..TIMER # The global module timer.
using ..Data
import ..build_model
import ..CutPattern, ..get_cut_pattern
import ..BuildStopReason, ..BUILT_MODEL, ..FOUND_OPTIMUM
using ..Utilities
using ..Utilities.Args
using ..PPG2KP, ..PPG2KP.Args
using ..Flow, ..Flow.Args
#using ..KnapsackPlates, ..KnapsackPlates.Args
import ..solution_value, ..to_pretty_str, ..simplify!, ..to_tikz_picture
# Used inside read_build_solve_and_print in tandem with generic-time-limit.
import ..throw_if_timeout_now, ..TimeoutError
# Used for SaveModel.write_to_file that, differently from JuMP.write_to_file,
# tries to use the raw solver to generate the MPS.
import ..SaveModel
const SM = SaveModel

"""
    create_unprefixed_subset(prefix, p_args :: T) :: T

!!! **Internal use.**

Given some `prefix`, query the `accepted_arg_list(Val{Symbol(prefix)})`,
to know which arguments were prefixed this way, search for them (with
the prefix) in `p_args` and return a new `typeof(p_args)` object in
which there is only the searched key-value pairs but the keys are
changed to not have the prefix anymore.
"""
function create_unprefixed_subset(prefix, p_args :: T) :: T where {T}
	subset = empty(p_args)
	prefix_id = Val(Symbol(prefix))
	prefix_options_names = getfield.(accepted_arg_list(prefix_id), :name)
	for name in prefix_options_names
		prefixed_name = prefix * "-" * name
		@assert haskey(p_args, prefixed_name)
		subset[name] = p_args[prefixed_name]
	end
	subset
end

"""
    round_instance(instance, factor, roundmode)

Return a new instance with all size-related fields multiplied by `factor`.

The `roundmode` (a `Base.RoundingMode`) specifies how the instances should
be rounded back to their original types (as the `factor` is expected to
be a `Float64` and not an `Integer`).

This function should be extended for new instance types; if you
want to use `--round-{nearest,up,down}` flags of
`GuillotineModels.CommandLine.run` for new instance types just create a
method that takes an instance of the new type, a `Float64`, and a
`Base.RoundingMode`).

This function is called with this signature only if at least one option
among `round-{nearest,up,down}` has a value different than one.

!!! Unmodified arrays (like `d`) may be shared between old and new instances.
"""
function round_instance(
	instance :: G2KP{D, S, P}, factor, roundmode
) where {D, S, P}
	@unpack L, W, l, w, d, p = instance
	L_, W_, l_, w_ = _round_instance(L, W, l, w, factor, roundmode)
	return G2KP{D, S, P}(L_, W_, l_, w_, d, p)
end

function round_instance(
	instance :: SLOPP{D, S, P}, factor, roundmode
) where {D, S, P}
	@unpack L, W, l, w, dlb, dub, p = instance
	L_, W_, l_, w_ = _round_instance(L, W, l, w, factor, roundmode)
	return SLOPP{D, S, P}(L_, W_, l_, w_, dlb, dub, p)
end

function round_instance(
	instance :: SSSCSP{D, S, P}, factor, roundmode
) where {D, S, P}
	@unpack L, W, l, w, d = instance
	L_, W_, l_, w_ = _round_instance(L, W, l, w, factor, roundmode)
	return SSSCSP{D, S, P}(L_, W_, l_, w_, d)
end

function _round_instance(
	L :: S, W :: S, l :: Vector{S}, w :: Vector{S}, factor, roundmode
) where {S}
	L_ = convert(S, round(L * factor, roundmode)) :: S
	W_ = convert(S, round(W * factor, roundmode)) :: S
	l_ = convert.(S, round.(l .* factor, roundmode)) :: Vector{S}
	w_ = convert.(S, round.(w .* factor, roundmode)) :: Vector{S}
	return L_, W_, l_, w_
end

"""
    round_instance(instance, p_args) -> old_or_new_instance

!!! **Internal use.**

Given a recognized instance type, uses `p_args` keys
`round-{nearest,up,down}` to either: (1) return them unmodified
if all keys have value one; (2) return a copy of them that is multiplied
by the ratio and rounded the specified way (no two keys may have a value
different than one).
"""
@timeit TIMER function round_instance(instance, p_args)
	# assert explanation: at least two of the three flags are disabled (i.e.,
	# have value one)
	@assert sum(isone.((
		p_args["round-nearest"],
		p_args["round-up"],
		p_args["round-down"]
	))) >= 2

	if p_args["round-nearest"] != 1
		factor = p_args["round-nearest"] :: Float64
		roundmode = RoundNearest
	elseif p_args["round-up"] != 1
		factor = p_args["round-up"] :: Float64
		roundmode = RoundUp
	elseif p_args["round-down"] != 1
		factor = p_args["round-down"] :: Float64
		roundmode = RoundDown
	else
		factor = 1
	end
	if factor != 1
		instance = round_instance(
			instance, factor, roundmode
		) :: typeof(instance)
	end
	return instance
end

# Fallback for unknown instance types.
function print_instance_stats(instance, verbose :: Bool)
	println("Start of the fallback output of `print_instance_stats`.")
	show(IOContext(stdout, :limit => false), instance)
	println("\nEnd of the fallback output of `print_instance_stats`.")
	return
end

function print_instance_stats(
	instance :: G2KP{D, S, P}, verbose :: Bool
) where {D, S, P}
	@unpack L, W, l, w, d = instance
	p_ = max(L/minimum(l), W/minimum(w))
	n_ = sum(d)
	n = length(d)
	if verbose
		@show L
		@show W
		@show n_
		@show n
		println("min_l = $(minimum(l))")
		println("min_w = $(minimum(w))")
		@show p_
		show(IOContext(stdout, :limit => false), instance)
		println()
	end
	return
end

# Get every $<something> in which the $ is not preceded by @
const TEMPLATE_REGEX = r"(?<!@)\$<([^>]+)>"
const ESCAPED_REGEX = r"@(\$<[^>]+>)"

# Treats `instance_file` specially.
function expand_template(
	template :: AbstractString, dict :: Dict{<:AbstractString, <:Any}
)
	expanded = replace(template, TEMPLATE_REGEX => function(s)
		# The below is byte indexing, not character indexing, but it is safe for
		# UFT-8 strings because '$<' are two bytes and '>' is another byte.
		key = s[3:end-1]
		if key == "instance_file"
			basename(dict["instance_path"])
		elseif haskey(dict, key)
			value = dict[key]
			if isa(value, Bool)
				string(convert(Int, value))
			else
				string(value)
			end
		else
			error("The template '\$<$key>' is invalid. There is no such parameter.")
		end
	end)
	cleaned = replace(expanded, ESCAPED_REGEX => s"\1")

	return cleaned
end

# Serves both as a function barrier and the step to be executed for each
# instance in the instance(s) file.
function build_solve_and_print(problem, formulation, instance_, pp, timings)
	start = time() :: Float64 # seconds since epoch
	limit = pp["generic-time-limit"] :: Float64
	verbose = !(pp["no-csv-output"] :: Bool)

	# We replace any $<parameter> pattern in `pp[save-model]` as soon as
	# possible, so if `parameter` is written wrong, the error is thrown early.
	if !isempty(pp["save-model"]) && verbose
		save_model_expanded = expand_template(pp["save-model"], pp)
		@show save_model_expanded
	end

	verbose && append!(timings, ["build_and_solve_time", "build_time"])
	instance = round_instance(instance_, pp)
	print_instance_stats(instance, verbose)

	throw_if_timeout_now(start, limit)

	solver_pp = create_unprefixed_subset(pp["solver"], pp)
	m = empty_configured_model(Val(Symbol(pp["solver"])), solver_pp)
	throw_if_timeout_now(start, limit)

	model_pp = create_unprefixed_subset(pp["model"], pp)
	bsr, mbp = build_model(
		problem, formulation, instance, m, model_pp
	)
	if verbose
		println("build_stop_reason = $(bsr)")
		println("build_stop_code = $(Int(bsr))")
	end

	if bsr == BUILT_MODEL
		if pp["relax2lp"]
			JuMP.relax_integrality(m)
		end
		num_vars = num_variables(m)
		num_constrs = num_all_constraints(m)
		if verbose
			@show num_vars
			@show num_constrs
		end
	end

	verbose && close_and_print!(timings, "build_time")
	throw_if_timeout_now(start, limit)

	if !isempty(pp["save-model"]) && verbose
		if bsr == BUILT_MODEL
			save_model_time = @elapsed begin
				if pp["use-native-save-model-if-possible"]
					@timeit TIMER "save_model" SM.write_to_file(m, save_model_expanded)
				else
					@timeit TIMER "save_model" JuMP.write_to_file(m, save_model_expanded)
				end
			end
			@show save_model_time
		else
			@warn "No model saved to '$(save_model_expanded)' because the" *
				" reason for stopping the model building was not $(BUILT_MODEL)" *
				" but instead $(bsr)."
		end
	end

	if pp["do-not-solve"]
		verbose && close_and_print!(timings, "build_and_solve_time")
		return nothing
	end

	throw_if_timeout_now(start, limit)

	if bsr == BUILT_MODEL
		verbose && println("MARK_FINAL_GENERIC_SOLVE")
		fms_name = "finished_model_solve"
		fms_time = @elapsed begin
			# The version taking an amount of seconds does not throw a timeout
			# error, it just sets the time limit.
			@timeit TIMER fms_name optimize_within_time_limit!(
				m, limit - (time() - start)
			)
		end
		if verbose
			println("$fms_name = $fms_time")
			stop_reason = termination_status(m)
			@show stop_reason
			stop_code = Int64(stop_reason)
			@show stop_code
		end
	end
	if verbose
		close_time = close_and_print!(timings, "build_and_solve_time")
		push!(timings, TimeSection("solution_print_time", close_time))
	end

	# This needs to be done even if it is not printed to warm-start the JIT
	# when a mock run is done.
	@assert bsr in (BUILT_MODEL, FOUND_OPTIMUM)
	solution = if bsr == BUILT_MODEL
		if primal_status(m) == MOI.FEASIBLE_POINT
			obj_value = objective_value(m)
			verbose && @show obj_value
			obj_bound = objective_bound(m)
			verbose && @show obj_bound
			# If the model was solved relaxed we try not to obtain a solution from
			# it (because it would not be a valid CutPattern). Also, we do not
			# print any info because we are formulation agnostic, and without
			# knowledge of the formulation, the value of the variables is of little
			# relevance.
			if pp["relax2lp"]
				nothing
			else
				get_cut_pattern(problem, formulation, m, mbp)
			end
		else
			nothing
		end
	elseif bsr == FOUND_OPTIMUM
		get_cut_pattern(problem, formulation, m, mbp)
	end
	if solution !== nothing
		@timeit TIMER "stringfy_solutions" begin
			solution_str = "solution = $solution\n"
			local pretty_sol_str :: String
			local pretty_simple_sol_str :: String
			if isa(solution, CutPattern)
				pretty_sol_str = to_pretty_str(solution)
				tikz_pretty_sol_str = to_tikz_picture(solution)

				simplified_sol = simplify!(deepcopy(solution))
				pretty_simple_sol_str = to_pretty_str(simplified_sol)
				tikz_pretty_simple_sol_str = to_tikz_picture(simplified_sol)
			else
				pretty_sol_str = join(to_pretty_str.(solution), "\n")
				tikz_pretty_sol_str = join(to_tikz_picture.(solution), "\n")

				simplified_sols = simplify!.(deepcopy(solution))
				pretty_simple_sol_str = join(to_pretty_str.(simplified_sols), "\n")
				tikz_pretty_simple_sol_str = join(
					to_tikz_picture.(simplified_sols), "\n"
				)
			end

			pretty_sol_str = "PRETTY_STR_SOLUTION_BEGIN\n" * pretty_sol_str *
				"\nPRETTY_STR_SOLUTION_END\n"
			tikz_pretty_sol_str = "TIKZ_PRETTY_STR_SOLUTION_BEGIN\n" *
				tikz_pretty_sol_str * "\nTIKZ_PRETTY_STR_SOLUTION_END\n"

			pretty_simple_sol_str = "SIMPLIFIED_PRETTY_STR_SOLUTION_BEGIN\n" *
				pretty_simple_sol_str *
				"\nSIMPLIFIED_PRETTY_STR_SOLUTION_END\n"
			tikz_pretty_simple_sol_str =
				"TIKZ_SIMPLIFIED_PRETTY_STR_SOLUTION_BEGIN\n" *
				tikz_pretty_simple_sol_str *
				"\nTIKZ_SIMPLIFIED_PRETTY_STR_SOLUTION_END\n"
		end

		if verbose
			@timeit TIMER "print_solutions" begin
				iob = IOBuffer()
				write(iob, solution_str)
				write(iob, pretty_sol_str)
				write(iob, tikz_pretty_sol_str)
				write(iob, pretty_simple_sol_str)
				write(iob, tikz_pretty_simple_sol_str)
				print(read(seekstart(iob), String))
			end
			sol_val = solution_value(problem, instance, solution)
			println("solution_value = $sol_val")
		end
	end
	verbose && close_and_print!(timings, "solution_print_time")

	# This is done just before returning because the version of
	# `optimize_within_time_limit!` called above does not throw. The
	# non-throwing version was used to be able to print as much information as
	# possible before throwing. However, if the timeout was reached, this
	# method MUST throw.
	if termination_status(m) == MOI.TIME_LIMIT
		throw(TimeoutError(start, limit, time()))
	end
	throw_if_timeout_now(start, limit)

	return
end

# Read the instance, build the model, solve the model, and print related stats.
"""
    read_build_solve_and_print(problem, format, pp)

!!! **Internal use.**

Given the parsed parameters (`pp`), read the instance file, build the model,
solve the model (unless `pp['do-not-solve']` is true), and print statistics
related to this process (unless `pp['no-csv-output']` is true).
The list of options recognized and implemented by this method is the
list returned by `generic_args()` (it also needs the required arguments
in `core_argparse_settings()`). The other arguments are solver or model
specific and are extracted and passed to their specific methods.
"""
@timeit TIMER function read_build_solve_and_print(
	format, pp
) # pp stands for parsed parameters
	verbose = !(pp["no-csv-output"] :: Bool)
	instance_path = pp["instance_path"] :: String
	verbose && @show instance_path
	timings = TimeSection[]
	verbose && push!(timings, "read_instance_time")
  dataset = Data.read_from_file(format, instance_path)
	verbose && close_and_print!(timings, "read_instance_time")
	problem = Val(Symbol(pp["problem"] :: String))
	model_type = Val(Symbol(pp["model"] :: String))

	if is_collection(format)
		for (index, instance) in enumerate(dataset)
			if verbose
				println("INSTANCE_START_MARKER_$index")
				push!(timings, "total_instance_time")
			end
			build_solve_and_print(problem, model_type, instance, pp, timings)
			if verbose
				close_and_print!(timings, "total_instance_time")
				println("INSTANCE_CLOSE_MARKER_$index")
			end
		end
	else
		if verbose
			println("INSTANCE_START_MARKER_1")
			verbose && push!(timings, "total_instance_time")
		end
		build_solve_and_print(problem, model_type, dataset, pp, timings)
		if verbose
			verbose && close_and_print!(timings, "total_instance_time")
			println("INSTANCE_CLOSE_MARKER_1")
		end
	end

	return nothing
end

"""
    gen_prefixed_argparse_settings(solver_or_model_name) :: ArgParseSettings

!!! **Internal use.**

Builds and returns an ArgParseSettings object representing all options
of the given solver or model already prefixed with its name.

The `solver_or_model_name` may be a Symbol or String (none is prefered),
and the method works because the solver or model implements
`accepted_arg_list(Val{Symbol(solver_or_model_name)})`
"""
function gen_prefixed_argparse_settings(
	solver_or_model_name :: Union{String, Symbol}
) :: ArgParse.ArgParseSettings
	prefix = solver_or_model_name # long name only for documentation
	s = ArgParseSettings()
	ArgParse.add_arg_group!(
		s, "$(prefix)-specific Options", "$(prefix)-specific-options"
	)
	original_args = accepted_arg_list(Val(Symbol(prefix)))
	prefixed_args = map(original_args) do arg
		Arg(string(prefix) * "-" * arg.name, arg.default, arg.help)
	end
	for arg in prefixed_args
		ArgParse.add_arg_table!(s, arg)
	end
	set_default_arg_group!(s)
	s
end

"""
    core_argparse_settings() :: ArgParseSettings

!!! **Internal use.**

An ArgParseSettings with the three core positional arguments `model`,
`solver`, and `instance_path`. They cannot be modeled as `Arg` objects because,
by design, all extra arguments must be options (i.e., be optional and preceded
by dashes).
"""
function core_argparse_settings() :: ArgParseSettings
	s = ArgParse.ArgParseSettings()
	ArgParse.add_arg_group!(s, "Core Parameters", "core_parameters")
	@add_arg_table! s begin
		"problem"
			help = "The type of problem to be solved (case sensitive, ex.: G2KP, G2CSP, G2ODP, G2MKP). Required."
			arg_type = String
		"format"
			help = "The format of the instance(s) in the file (case sensitive, ex.: Classic_G2KP, CPG_SLOPP, CPG_SSSCSP, CPG_ODPW, CPG_MHLOPPW). Required."
			arg_type = String
		"model"
			help = "Model or solution procedure to be used (case sensitive, ex.: Flow, KnapsackPlates, PPG2KP). Required."
			arg_type = String
		"solver"
			help = "Solver to be used if necessary (case-sensitive, ex.: NoSolver, Cbc, CPLEX, Gurobi). Required, even if --do-not-solve is specified. NoSolver use `JuMP.Model()` which allows building and saving the model, but not solving it."
			arg_type = String
		"instance_path"
			help = "The path to the instance to be solved."
			arg_type = String
	end
	set_default_arg_group!(s)
	s
end

"""
    generic_args() :: Vector{Arg}

!!! **Internal use.**

All the `Arg`s representing options that are independent from solver or model.
"""
function generic_args() :: Vector{Arg}
	return [
		Arg(
			"do-not-solve", false,
			"The model is built but not solved. A solver has yet to be specified. Note that just building a model may depend on using a solver over subproblems. Such uses of the solver are not disabled by this flag."
		),
		Arg(
			"generic-time-limit", 60.0 * 60.0 * 24.0 * 365.0,
			"Defines a time limit (in seconds) to be observed in the context of the model-agnostic process of reading, building, solving, and printing. Each model has to define its own flag to control the time inside the model building process (to be set independently, it is does not interact with this flag). `JuMP.set_time_limit_sec` is called over the model before starting to solve it, with the remaining time after reading instance and building the model (not the total time). If `--warm-jit` defines that there will be a mock run to warm the jit, the time limit applies to this mock run (i.e., if the mock run breaks the limit it an exception is raised) but the timer is reset after the mock, so the mock time is not counted for the 'real' run time limit. A `GuillotineModels.TimeoutError` is raised if the time limit is violated."
		),
		Arg(
			"save-model", "",
			"Save the model of the passed instance to the given filename. The format used depends on the extension of the filename, allowed formats are described by the enumeration `MathOptInterface.FileFormats.FileFormat` from package `MathOptInterface`). The filename also allows using \$<parameter_name> patterns inside the name, these will be replaced by the value of the parameter (if it was not passed, the default value will be used). For convenience, a pseudo-parameter instance_file is also available (it is the same as `basename(instance_path)`, it includes any file extensions), and `Bool` parameters (i.e., flags with no argument) are replaced by 0 or 1. Putting an @ in front of \$<some_string> will disable the substitution (and remove the @ from the final string). There is no way to express an @ followed by a substitution that actually occurs (this is a limitation of the code)."
		),
		Arg(
			"use-native-save-model-if-possible", false,
			"Only used if a filename is passed to '--save-model'. If `--use-native-save-model-if-possible` is passed, then `GuillotineModels.SaveModel.write_to_file` is used instead of `JuMP.write_to_file`. Check the `SaveModel` module documentation for more info. For at least Gurobi and CPLEX, `GuillotineModels.SaveModel.write_to_file` will use the raw solver model-saving utilities instead of the `JuMP` one, but if the solver is not recognized then it will default to the JuMP one."
		),
		Arg(
			"warm-jit", "no",
			"There are three valid options: 'no', 'with-toy', and 'yes'. If 'yes', the instance is solved two times, but the first time the arguments are changed to disable output (i.e., '--no-csv-output' and '--SOLVER-no-output' are passed, and if supported, '--MODEL-quiet' is passed too; in the GuillotineModels.TIMER the timings of this first run are under 'warm_jit'). If 'with-toy', it behaves similarly to 'yes' but instead of using the instance itself it uses a very small hardcoded instance (this helps a lot, but many procedures only called when the model has specific properties are not called, and therefore, not compiled). If 'no', the code is just run a single time."
		),
		Arg(
			"no-csv-output", false,
			"Disable some extra output that the author of this package uses to assemble CSVs and/or debug. Also disables --save-model."
		),
		Arg(
			"relax2lp", false,
			"Integer and binary variables become continuous."
		),
		Arg(
			"round-nearest", 1.0,
			"Multiplies the instances size-related fields by the passed factor and round them to nearest (the model answer becomes a GUESS, not a valid primal heuristic, nor a valid bound)."
		),
		Arg(
			"round-up", 1.0,
			"Multiplies the instances size-related fields by the passed factor and round them up (the model becomes a PRIMAL HEURISTIC)."
		),
		Arg(
			"round-down", 1.0,
			"Multiplies the instance size-related fields by the passed factor and round them down (the model becomes an OPTIMISTIC GUESS, A VALID BOUND)."
		)
	]
end

"""
    generic_argparse_settings() :: ArgParseSettings

!!! **Internal use.**

An ArgParseSettings containing all options (not positional arguments) that are
independent from the chosen solver or model.
"""
function generic_argparse_settings()
	s = ArgParse.ArgParseSettings()
	add_arg_group!(s, "Generic Options", "generic_options")
	for arg in generic_args()
		add_arg_table!(s, arg)
	end
	set_default_arg_group!(s)
	s
end

"""
    argparse_settings(models_list, solvers_list) :: ArgParseSettings

!!! **Internal use.**

Create an ArgParseSettings which includes the core arguments, generic
arguments, and all arguments from models and solvers available (prefixed by
their model or solver name).

Options of all models and solvers are added even if just one model is selected
at a time. The reasons for this are: (1) we do not know the model before the
parsing unless we manipulate ARGS directly; (2) if the unused solvers and
models are not included in the ArgParseSettings they do not appear in the help
message.
"""
function argparse_settings(
	models_list :: Vector{Symbol}, solvers_list :: Vector{Symbol}
) :: ArgParseSettings
	s = core_argparse_settings()
	ArgParse.import_settings!(s, generic_argparse_settings())
	for sym in [solvers_list; models_list]
		ArgParse.import_settings!(s, gen_prefixed_argparse_settings(sym))
	end
	s
end

"""
    warn_if_changed_unused_values(p_args, models_list, solvers_list)

!!! **Internal use.**

Gives warning messages if `p_args` has a value different from the default
for an option of a model or solver that is not the one used. In other words,
help a distracted user to not keep thinking it has passed an option to the
used solver/model when they have not. Many solvers have the same option
but with a different prefix, it is easy to change the solver used and forget
to also change the prefix in the parameter options.

It is actually impossible to know if an argument was passed or not by the
command-line after they are parsed by ArgParse because all of the arguments
must have a default value and if the argument is not present in the
command-line their default value is placed in the dict by ArgParse. This is
the reason we just check if non-used options have values different from
default, instead of just checking if they were provided.

This should work out-of-the-box for any third-party models or solvers
given they implement their own version of `accepted_arg_list` and also
have their identifying symbol passed in either `models_list` or
`solvers_list`.
"""
function warn_if_changed_unused_values(p_args, models_list, solvers_list)
	# Remember that for now, only solver and model support options/flags.
	# problem and format do not support them yet.
	core_arg_names = ["problem", "format", "solver", "model", "instance_path"]
	generic_arg_names = getfield.(generic_args(), :name)

	model = p_args["model"]
	model_id = Val(Symbol(model))
	model_arg_names = getfield.(accepted_arg_list(model_id), :name)
	model_arg_names = map(name -> model * "-" * name, model_arg_names)
	# the same as the block above just with "model" replaced by "solver"
	solver = p_args["solver"]
	solver_id = Val(Symbol(solver))
	solver_arg_names = getfield.(accepted_arg_list(solver_id), :name)
	solver_arg_names = map(name -> solver * "-" * name, solver_arg_names)

	expected_arg_names = [
		model_arg_names; core_arg_names; generic_arg_names; solver_arg_names
	]
	sort!(expected_arg_names)

	used_groups = Symbol.([p_args["model"], p_args["solver"]])
	unused_groups = setdiff([models_list; solvers_list], used_groups)
	unused_prefixed_args = Arg[]
	for unused_group in unused_groups
		args = accepted_arg_list(Val(unused_group))
		prefixed_args = map(args) do arg
			Arg(string(unused_group) * "-" * arg.name, arg.default, arg.help)
		end
		append!(unused_prefixed_args, prefixed_args)
	end
	unused_prefixed_args = collect(Iterators.flatten(unused_prefixed_args))

	for arg_name in keys(p_args)
		if isempty(searchsorted(expected_arg_names, arg_name))
			unused_arg_idx = findfirst(a -> a.name == arg_name, unused_prefixed_args)
			if unused_arg_idx === nothing
				@error(
					"ArgParse accepted option --$(arg_name) but it is not a" *
					" recognized option. This should not be possible."
				)
			else
				unused_arg = unused_prefixed_args[unused_arg_idx]
				unused_arg.default == p_args[arg_name] && continue
				@warn(
					"Option --$(unused_arg.name) has a different value than the" *
					" default, but it is not from the solver or model selected." *
					" You probably did use an option of an unselected solver or model."
				)
			end
		end
	end

	return nothing
end

"""
    throw_if_incompatible_options(p_args)

!!! **Internal use.**

Check the already parsed arguments and test if options that are incompatible
with each other were used.
"""
function throw_if_incompatible_options(p_args)
	# Generic Flags conflicts
	num_rounds = sum(.! isone.(map(flag -> p_args[flag],
		["round-nearest", "round-up", "round-down"]
	)))
	num_rounds > 1 && @error(
		"only one of --round-{nearest,up,down} may be passed at the" *
		" same time (what the fuck you expected?)"
	)
	# The flag conflicts of the specific model and solver.
	Utilities.Args.throw_if_incompatible_options(
		Val(Symbol(p_args["model"])),
		create_unprefixed_subset(p_args["model"], p_args)
	)
	Utilities.Args.throw_if_incompatible_options(
		Val(Symbol(p_args["solver"])),
		create_unprefixed_subset(p_args["solver"], p_args)
	)
	throw_if_unrecognized(
		"warm-jit", p_args["warm-jit"], ("no", "with-toy", "yes")
	)
	option_name = "generic-time-limit"
	p_args[option_name] > 0.0 || throw(ArgumentError(
		"Option $(option_name) must be positive, but it is" *
		" $(p_args[option_name])."
	))
end

# Definition of the command line arguments.
"""
    parse_args(args, models_list, solvers_list) :: Dict{String, Any}

!!! **Internal use.**

Given a vector of the command-line arguments `args` and the lists of
available models and solvers, parse the arguments. If the `args`
refer to a model or solver not in `models_list` or `solvers_list`
exceptions may be thrown. If `args` just triggers the help message,
an empty `Dict` is returned.
"""
@timeit TIMER function parse_args(
	args, models_list, solvers_list
) :: Dict{String, Any}
	s = argparse_settings(models_list, solvers_list)
	p_args = ArgParse.parse_args(args, s)
	isnothing(p_args) && return Dict{String, Any}()
	warn_if_changed_unused_values(p_args, models_list, solvers_list)
	throw_if_incompatible_options(p_args)

	return p_args
end

"""
    toy_instance(problem, format):: String

Return a `String` representing an instance of problem in `format`.

This is the function that should be extended in order to be able to call `run`
with the `--warm-jit with-toy`.
"""
function toy_instance(problem, format)
	error(
		"If you want to mock a run with `--warm-jit with-toy` there need to" *
		" exist an implementation of `GuillotineModels.CommandLine." *
		"toy_instance` taking the instance format you passed. You have passed" *
		" problem `$problem` and format `$format`."
	)
end

function toy_instance(::Val{:G2KP}, ::Val{:CPG_SLOPP}) :: String
	return toy_instance(Val(:G2KP), CPG_SLOPP{Int,Int,Int}())
end

function toy_instance(
	::Val{:G2KP}, ::CPG_SLOPP{D, S, P}
) :: String where {D, S, P}
"""
***2D Rectangular Problem***
***Instances for the Single Large Object Placement Problem (SLOPP)***
Input parameter file: SLOPP_parameters.txt
***************************************************************************************************************
Total number of instances
LargeObject.Length      LargeObject.Width
Number of different item types (i)
Item[i].Length  Item[i].Width   Item[i].LowerBoundDemand        Item[i].UpperBoundDemand        Item[i].Value
***************************************************************************************************************
1
100	200
2
50     100    0     1     2
50     200    0     1     1
"""
end

function toy_instance(::Val{:G2MKP}, ::Val{:CPG_MHLOPPW}) :: String
	return toy_instance(Val(:G2MKP), CPG_MHLOPPW{Int, Int, Int}())
end
function toy_instance(
	::Val{:G2MKP}, ::CPG_MHLOPPW{D, S, P}
) :: String where {D, S, P}
"""
***2D Rectangular Problem***
***Instances for the Multiple Heterogeneous Large Object Placement Problem (MHLOPP/W)***
Input parameter file: written_by_hand
***************************************************************************************************************
Total number of instances
Number of different large objects (j)
LargeObject[j].Length	LargeObject[j].Width	LargeObject[j].Available
Number of different item types (i)
Item[i].Length	Item[i].Width	Item[i].LowerBoundDemand	Item[i].UpperBoundDemand	Item[i].Value
***************************************************************************************************************
1
1
100	200	2
4
100	200	0	3	20000
50	150	0	1	7501
50	200	0	2	10001
50	50	0	2	2501"""
end

const CPG_SSSCSP_HEADER = """
***2D Rectangular Problem***
***Instances for the Single Stock Size Cutting Stock Problem (SSSCSP)***
Input parameter file: written_by_hand
****************************************************************************************************
Total number of instances
LargeObject.Length	LargeObject.Width
Number of different item types (i)
Item[i].Length	Item[i].Width	Item[i].Demand
*****************************************************************************************************"""

function toy_instance(::Val{:G2OPP}, ::Val{:CPG_SSSCSP}) :: String
	return toy_instance(Val(:G2OPP), CPG_SSSCSP{Int, Int, Int}())
end
function toy_instance(
	::Val{:G2OPP}, ::CPG_SSSCSP{D, S, P}
) :: String where {D, S, P}
"""
$CPG_SSSCSP_HEADER
1
100	200
3
100	50	2
25	100	2
50	100	1"""
end

# Same as G2OPP, but with the double of the items and therefore needing
# two original plates (i.e., large objects).
function toy_instance(::Val{:G2CSP}, ::Val{:CPG_SSSCSP}) :: String
	return toy_instance(Val(:G2CSP), CPG_SSSCSP{Int, Int, Int}())
end
function toy_instance(
	::Val{:G2CSP}, ::CPG_SSSCSP{D, S, P}
) :: String where {D, S, P}
"""
$CPG_SSSCSP_HEADER
1
100	200
3
100	50	4
25	100	4
50	100	2"""
end

function toy_instance(::Val{:G2KP}, ::Val{:Classic_G2KP}) :: String
	return toy_instance(Val(:G2KP), Classic_G2KP{Int,Int,Int}())
end

function toy_instance(
	::Val{:G2KP}, ::Classic_G2KP{D, S, P}
) :: String where {D, S, P}
"""
100	200
2
50     100    1     2
50     200    1     1
"""
end

"""
    warm_jit(p_args)

!!! **Internal use.**

Looks at `p_args["warm-jit"]` and implements what is said in the flag
description.

See also: [`toy_instance`](@ref), [`read_build_solve_and_print`](@ref)
"""
@timeit TIMER function warm_jit(p_args)
	p_args["warm-jit"] == "no" && return
	problem = Val(Symbol(p_args["problem"] :: String))
	format = Val(Symbol(p_args["format"] :: String))
	copyed_args = deepcopy(p_args)
	# The temporary file is only needed if p_args["warm-jit"] == "with-toy",
	# however, if it is needed, then it needs to be valid for the rest of
	# this method scope.
	mktemp() do path, io
		if p_args["warm-jit"] == "with-toy"
				write(io, toy_instance(problem, format) :: String)
				close(io)
				copyed_args["instance_path"] = path
		end
		copyed_args["no-csv-output"] = true
		copyed_args[p_args["solver"] * "-" * "no-output"] = true
		valid_model_args = accepted_arg_list(Val(Symbol(p_args["model"])))
		if any(arg -> arg.name == "quiet", valid_model_args)
			copyed_args[p_args["model"] * "-" * "quiet"] = true
		end
		read_build_solve_and_print(format, copyed_args)
	end
	return
end

# TODO: For now, there is no help support for which are possible values of
# the problem and the format positional parameters. There is also no way
# to support problem-specific (or format-specific) options/flags.
"""
    run(args = ARGS; implemented_models = [...], supported_solvers = [...])

Parse the command-line arguments and take appropriate action: prints help
or solve instance(s) in the specified file of the specified format for
the specified problem using the specified solver, model, and their options.

The parameters available (listed in the help message and actually parsed)
depend on the `implemented_models` and `supported_solvers`. The default values
are the models and solvers made available by the GuillotineModels package. If
you implement your own models or add support for more solvers you need to pass
them there for them to be considered (and if you want the old ones to keep
working you need to specify them also).

The best way to know everything this command is capable is to call:
```
import SUPPORTED_SOLVERS_YOU_HAVE_AVAILABLE
import GuillotineModels
GuillotineModels.CommandLine.run(
	["--help"];
	supported_solvers = [SUPPORTED_SOLVERS_YOU_HAVE_AVAILABLE]
)
```

Which will give you the help message. If you want to use it as an script
you just need to remove the `["--help"]` from the call.
"""
@timeit TIMER function run(
	args = ARGS;
	implemented_models :: Vector{Symbol} = [:Flow, :PPG2KP],#, :KnapsackPlates],
	supported_solvers :: Vector{Symbol} = [
		:CPLEX, :Gurobi, :Cbc, :GLPK, :NoSolver
	]
)
	p_args = parse_args(args, implemented_models, supported_solvers)
	isempty(p_args) && return # Happens if called just for "--help".
	verbose = !(p_args["no-csv-output"] :: Bool)
	verbose && @show args
	verbose && @show p_args
	# Out of the `if verbose` because we want it to be compiled in mock runs.
	date_now = Dates.format(Dates.now(), dateformat"yyyy-mm-ddTHH:MM:SS")
	verbose && @show date_now
	warm_jit(p_args) # It checks `p_args["warm-jit"]` to check if it will run.
	# remaining code should not depend on this option
	delete!(p_args, "warm-jit")
	# Now that is safe to extract some options, let us already call the
	# read_build_solve_and_print function with the correct format,
	# so it can be specialized for it.
	format = Val(Symbol(p_args["format"] :: String))
	total_file_time = @elapsed read_build_solve_and_print(format, p_args)
	verbose && @show total_file_time

	return
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 59aaecb0-bde0-11eb-3d1a-496b1d550fba
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="DifferentialEquations", version="6"),
	    Pkg.PackageSpec(name="Parameters", version="0.12"),
        Pkg.PackageSpec(name="LabelledArrays", version="1"),
        Pkg.PackageSpec(name="Setfield", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra, DifferentialEquations, Parameters, LabelledArrays, Setfield

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ e04166d2-a0b8-4e08-ad95-f16a151b0aba
md"""
# Figure 3.03 Michaelis-Menten kinetics
"""

# ╔═╡ 8c64a434-f663-4a26-85d3-326c6e57b1dc
# Full enzyme catalysis model
function full_model!(du, u, p, t)
	@unpack ET, K1, KM1, K2 = p
	@unpack S, ES, P = u
	freeEnz = ET - ES
	v1 = K1 * S * freeEnz - KM1 * ES
	v2 = K2 * ES
	du.S = - v1
	du.ES = v1 - v2
	du.P = v2
	return du
end

# ╔═╡ 594ae902-cc2a-42bd-b7c2-7852f66c49d1
u0 = LVector(S= 5.0, ES = 0.0, P = 0.0)

# ╔═╡ 7a2c5f66-b383-4ac8-9c09-425e9aeb0d2e
tend = 1.0

# ╔═╡ 52fde19e-0ed3-4472-81c0-37c4caca4afc
p = (ET = 1.0, K1 = 30.0, KM1 = 1.0, K2 = 10.0)

# ╔═╡ fc25f239-3f0d-4ca3-a356-95a282630f6c
sol = solve(ODEProblem(full_model!, u0, tend, p));

# ╔═╡ 47491dcd-7751-497e-bb9b-cb9097bae71c
begin
	p1 = plot(sol)
	plot!(p1, sol, vars=((t, es) -> (t, p.ET - es), 0, 2), label="E")
	plot!(p1, xlabel="Time (arbitrary units)", ylabel="Concentration (arbitrary units)", legend=:right)

	p1
end

# ╔═╡ ba148888-dc13-45e8-87a8-ce7bd9ecc93f
u0R = LVector(S = sum(u0))

# ╔═╡ 6056d577-a2ad-45d3-bd30-0b1df2eddbc9
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ f2ba7dc5-6556-4e4e-941b-7cc9de495e43
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
end

# ╔═╡ eb96706d-6f7f-4f35-b3ac-c5bfe0a71caf
# Apply QSSA on the ES complex to reduce the model complexity
function reduced_model!(du, u, p, t)
	@unpack ET, K1, KM1, K2 = p
	@unpack S = u
	du.S = -K2 * ET * hill(S, (KM1 + K2) / K1)
	return du
end

# ╔═╡ 884dd7d7-0ef7-4fb1-8208-74d7f4ded49d
solr = solve(ODEProblem(reduced_model!, u0R, tend, p));

# ╔═╡ 57a78828-f389-4b12-ae64-5c114e01401e
begin
	p2 = plot(sol, vars=(0, [1, 3]), line=(:dash), label=["S (full)", "P (full)"])
	plot!(p2, solr, line=(:blue), lab="S (reduced)")
	plot!(p2, solr, vars=((t, s)->(t, u0R[1] - s), 0, 1), line=(:red), lab="P (reduced)")
	plot!(p2, xlabel="Time (arbitrary units)",  ylabel="Concentration (arbitrary units)", xlims=(0.0,1.0), ylims=(0.0,5.0), legend = :right)

	p2
end

# ╔═╡ b4121e60-6726-42c1-a9fb-489bf988bc2d

# ╔═╡ Cell order:
# ╠═e04166d2-a0b8-4e08-ad95-f16a151b0aba
# ╠═8c64a434-f663-4a26-85d3-326c6e57b1dc
# ╠═eb96706d-6f7f-4f35-b3ac-c5bfe0a71caf
# ╠═594ae902-cc2a-42bd-b7c2-7852f66c49d1
# ╠═7a2c5f66-b383-4ac8-9c09-425e9aeb0d2e
# ╠═52fde19e-0ed3-4472-81c0-37c4caca4afc
# ╠═fc25f239-3f0d-4ca3-a356-95a282630f6c
# ╠═47491dcd-7751-497e-bb9b-cb9097bae71c
# ╠═ba148888-dc13-45e8-87a8-ce7bd9ecc93f
# ╠═884dd7d7-0ef7-4fb1-8208-74d7f4ded49d
# ╠═57a78828-f389-4b12-ae64-5c114e01401e
# ╠═6056d577-a2ad-45d3-bd30-0b1df2eddbc9
# ╠═f2ba7dc5-6556-4e4e-941b-7cc9de495e43
# ╠═59aaecb0-bde0-11eb-3d1a-496b1d550fba
# ╠═b4121e60-6726-42c1-a9fb-489bf988bc2d

ext: jl

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
# Exponential distribution

export Exponential

@parameterized Exponential(β) ≃ Lebesgue(ℝ₊)

@kwstruct Exponential()

function logdensity(d::Exponential{()} , x)
    return -x
end

Base.rand(rng::AbstractRNG, T::Type, μ::Exponential{()}) = randexp(rng,T)

TV.as(::Exponential) = asℝ₊

##########################
# Scale β

@kwstruct Exponential(β)

function Base.rand(rng::AbstractRNG, T::Type, d::Exponential{(:β,)})
    randexp(rng, T) * d.β
end

function logdensity(d::Exponential{(:β,)}, x)
    z = x / d.β
    return logdensity(Exponential(), z) - log(d.β)
end

distproxy(d::Exponential{(:β,)}) = Dists.Exponential(d.β)

asparams(::Type{<:Exponential}, ::Val{:β}) = asℝ₊

##########################
# Log-Scale logβ

@kwstruct Exponential(logβ)

function Base.rand(rng::AbstractRNG, T::Type, d::Exponential{(:logβ,)})
    randexp(rng, T) * exp(d.logβ)
end

function logdensity(d::Exponential{(:logβ,)}, x)
    z = x * exp(-d.logβ)
    return logdensity(Exponential(), z) - d.logβ
end

distproxy(d::Exponential{(:logβ,)}) = Dists.Exponential(exp(d.logβ))

asparams(::Type{<:Exponential}, ::Val{:logβ}) = asℝ

##########################
# Rate λ

@kwstruct Exponential(λ)

function Base.rand(rng::AbstractRNG, T::Type, d::Exponential{(:λ,)})
    randexp(rng, T) / d.λ
end

function logdensity(d::Exponential{(:λ,)}, x)
    z = x * d.λ
    return logdensity(Exponential(), z) + log(d.λ)
end

distproxy(d::Exponential{(:λ,)}) = Dists.Exponential(1/d.λ)

asparams(::Type{<:Exponential}, ::Val{:λ}) = asℝ₊

##########################
# Log-Rate logλ

@kwstruct Exponential(logλ)

function Base.rand(rng::AbstractRNG, T::Type, d::Exponential{(:logλ,)})
    randexp(rng, T) * exp(-d.logλ)
end

function logdensity(d::Exponential{(:logλ,)}, x)
    z = x * exp(d.logλ)
    return logdensity(Exponential(), z) + d.logλ
end

distproxy(d::Exponential{(:logλ,)}) = Dists.Exponential(exp(-d.logλ))

asparams(::Type{<:Exponential}, ::Val{:logλ}) = asℝ


text:
module Julia

include("pix.jl")
include("ppm.jl")

function julia(w, h)
  pix = Pix.zero(w,h)

  xa = -2.0
  xb = 2.0
  ya = -1.5
  yb = 1.5
  max_it = 65535

  # Find a good julia spot from the mandelbrot set.
  c = 0
  while true
    cx = rand() * (xb - xa) + xa
    cy = rand() * (yb - ya) + ya
    c = cx + cy * 1im
    z = c
    i = 0
    for i=0:max_it
      if abs(z) > 2.0
        break
      end
      z = z * z + c
    end
    if i > 10 && i < 32000
      break
    end
  end

  # Write the pixels
  for ih = 1:h
    zy = ih * (yb - ya) / (h - 1) + ya
    for iw = 1:w
      zx = iw * (xb - xa) / (w - 1) + xa
      z = zx + zy * 1im
      i = 0
      for i = 0:max_it
        if abs(z) > 2.0
          break
        end
        z = z * z + c
      end
      pix[iw,ih][1] = i % 128 * 512
      pix[iw,ih][2] = i % 256 * 256
      pix[iw,ih][3] = i % 512 * 128
    end
  end

  return pix
end

PPM.write(julia(1024,1024))

end

max_stars_repo_name: dominichamon/juliajulia

max_issues_repo_name: dominichamon/juliajulia

max_forks_repo_name: dominichamon/juliajulia


text:
struct Lengths{T}
    a::T # height of pivot
    b::T # length of long arm
    c::T # length of short arm
    d::T # lenght of weight arm
    e::T # length of sling
    u::T # radius of weight
    z::T # radius of projectile
end

struct Masses{T}
    w::T # mass of weight
    p::T # mass of projectile
    a::T # mass of arm
end

struct Angles{T} # in radians
    aq::T # angle between pivot stand and arm
    wq::T # angle between arm and weight
    sq::T # angle between arm and sling
end

struct AnglularVelocities{T}
    aw::T # for aq
    ww::T # for wq
    sw::T # for sq
end

struct Inertias{T}
    iw::T # inertia of weight
    ia::T # inertia of arm
end

struct Vec
    x
    y
end

struct Constants{T}
    w::T # wind speed
    ρ::T # Density of Air
    Cd::T # Drag Co-efficient of Projectile
    Grav::T
    r::T # release angle of projectile [input]
end

mutable struct Solution
    WeightCG
    WeightArm
    ArmSling
    Projectile
    SlingEnd
    ArmCG
    Time
    ReleaseVelocity
    ReleasePositon
end

mutable struct TrebuchetState
    l::Lengths
    m::Masses
    a::Angles # [stage 1 & 2]
    aw::AnglularVelocities # [stage 1 & 2]
    c::Constants
    i::Inertias
    stage::Union{Val{:Ground},Val{:Hang},Val{:Released},Val{:End}}
    rate
    p # projectile point [stage 3]
    v # projectile speed [stage 3]
    sol::Solution
    Tn # tension on string
    terminated
    function TrebuchetState(l::Lengths, m::Masses, c::Constants, rate)
        θ = asin(l.a/l.b)
        sq = π - θ
        aq = π/2 + θ
        wq = -aq
        ai = m.a*((l.b + l.c)^2)/12
        wi = m.w*(l.u)^2/2 # inertia for disk
        a = Angles(aq, wq, sq)
        aw = AnglularVelocities(0.0, 0.0, 0.0)
        i = Inertias(wi, ai)
        new(l, m, a, aw, c, i, Val{:Ground}(), rate, -1, -1, Solution(), 0, false)
    end
end


text:
module TestHash
using Test
using CategoricalArrays

@testset "hash() for CategoricalPool{Int} and CategoricalPool{Float64} and its values" begin
    pool1 = CategoricalPool([1, 2, 3])
    pool2 = CategoricalPool([2.0, 1.0, 3.0])

    @test (hash(pool1) == hash(pool1)) === true
    @test (hash(pool1) == hash(pool2)) === false
    @test (hash(pool2) == hash(pool2)) === true

    opool1 = CategoricalPool([1, 2, 3], true)
    opool2 = CategoricalPool([2.0, 1.0, 3.0], true)

    @test (hash(opool1) == hash(opool1)) === true
    @test (hash(opool1) == hash(opool2)) === false
    @test (hash(opool2) == hash(opool2)) === true

    nv1a = CategoricalValue(1, pool1)
    nv2a = CategoricalValue(1, pool2)
    nv1b = CategoricalValue(2, pool1)
    nv2b = CategoricalValue(2, pool2)

    @test (hash(nv1a) == hash(nv1a)) === true
    @test (hash(nv1a) == hash(nv2a)) === false
    @test (hash(nv1a) == hash(nv1b)) === false
    @test (hash(nv1a) == hash(nv2b)) === true

    @test (hash(nv1b) == hash(nv1a)) === false
    @test (hash(nv1b) == hash(nv2a)) === true
    @test (hash(nv1b) == hash(nv1b)) === true
    @test (hash(nv1b) == hash(nv2b)) === false

    @test (hash(nv2a) == hash(nv1a)) === false
    @test (hash(nv2a) == hash(nv2a)) === true
    @test (hash(nv2a) == hash(nv1b)) === true
    @test (hash(nv2a) == hash(nv2b)) === false

    @test (hash(nv2b) == hash(nv1a)) === true
    @test (hash(nv2b) == hash(nv2a)) === false
    @test (hash(nv2b) == hash(nv1b)) === false
    @test (hash(nv2b) == hash(nv2b)) === true

    ov1a = CategoricalValue(1, opool1)
    ov2a = CategoricalValue(1, opool2)
    ov1b = CategoricalValue(2, opool1)
    ov2b = CategoricalValue(2, opool2)

    @test (hash(ov1a) == hash(ov1a)) === true
    @test (hash(ov1a) == hash(ov2a)) === false
    @test (hash(ov1a) == hash(ov1b)) === false
    @test (hash(ov1a) == hash(ov2b)) === true

    @test (hash(ov1b) == hash(ov1a)) === false
    @test (hash(ov1b) == hash(ov2a)) === true
    @test (hash(ov1b) == hash(ov1b)) === true
    @test (hash(ov1b) == hash(ov2b)) === false

    @test (hash(ov2a) == hash(ov1a)) === false
    @test (hash(ov2a) == hash(ov2a)) === true
    @test (hash(ov2a) == hash(ov1b)) === true
    @test (hash(ov2a) == hash(ov2b)) === false

    @test (hash(ov2b) == hash(ov1a)) === true
    @test (hash(ov2b) == hash(ov2a)) === false
    @test (hash(ov2b) == hash(ov1b)) === false
    @test (hash(ov2b) == hash(ov2b)) === true

    @testset "ordered and non-ordered values hash equal" begin
    @test (hash(ov1a) == hash(nv1a)) === true
    @test (hash(ov1a) == hash(nv2a)) === false
    @test (hash(ov1a) == hash(nv1b)) === false
    @test (hash(ov1a) == hash(nv2b)) === true

    @test (hash(ov1b) == hash(nv1a)) === false
    @test (hash(ov1b) == hash(nv2a)) === true
    @test (hash(ov1b) == hash(nv1b)) === true
    @test (hash(ov1b) == hash(nv2b)) === false

    @test (hash(ov2a) == hash(nv1a)) === false
    @test (hash(ov2a) == hash(nv2a)) === true
    @test (hash(ov2a) == hash(nv1b)) === true
    @test (hash(ov2a) == hash(nv2b)) === false

    @test (hash(ov2b) == hash(nv1a)) === true
    @test (hash(ov2b) == hash(nv2a)) === false
    @test (hash(ov2b) == hash(nv1b)) === false
    @test (hash(ov2b) == hash(nv2b)) === true
    end
end

end


text:
function apply_prec_ilu0(bilu0, v)
  pv  = zeros( size(v)[1] )
  ccall( (:apply_prec_ILU0, "../LF/extlibs/libmkl.so"), Void,
         (Int, Ptr{Float64}, Ptr{Int}, Ptr{Int}, Ptr{Float64}, Ptr{Float64}),
          bilu0.n, bilu0.nzval, bilu0.colptr, bilu0.rowval, v, pv )
  return pv
end


text:
"""
    DiscreteUniform(a,b)

A *Discrete uniform distribution* is a uniform distribution over a consecutive sequence of integers between `a` and `b`, inclusive.

```math
P(X = k) = 1 / (b - a + 1) \\quad \\text{for } k = a, a+1, \\ldots, b.
```

```julia
DiscreteUniform(a, b)   # a uniform distribution over {a, a+1, ..., b}

params(d)       # Get the parameters, i.e. (a, b)
span(d)         # Get the span of the support, i.e. (b - a + 1)
probval(d)      # Get the probability value, i.e. 1 / (b - a + 1)
minimum(d)      # Return a
maximum(d)      # Return b
```

External links

* [Discrete uniform distribution on Wikipedia](http://en.wikipedia.org/wiki/Uniform_distribution_(discrete))
"""
struct DiscreteUniform <: DiscreteUnivariateDistribution
    a::Int
    b::Int
    pv::Float64 # individual probabilities

    function DiscreteUniform(a::Real, b::Real; check_args::Bool=true)
        check_args && @check_args(DiscreteUniform, a <= b)
        new(a, b, 1 / (b - a + 1))
    end
    DiscreteUniform(b::Real; check_args::Bool=true) = DiscreteUniform(0, b; check_args=check_args)
    DiscreteUniform() = new(0, 1, 0.5)
end

@distr_support DiscreteUniform d.a d.b

### Parameters

span(d::DiscreteUniform) = d.b - d.a + 1
probval(d::DiscreteUniform) = d.pv
params(d::DiscreteUniform) = (d.a, d.b)

### Show

show(io::IO, d::DiscreteUniform) = show(io, d, (:a, :b))

### Statistics

mean(d::DiscreteUniform) = middle(d.a, d.b)

median(d::DiscreteUniform) = fld(d.a + d.b, 2)

var(d::DiscreteUniform) = (span(d)^2 - 1.0) / 12.0

skewness(d::DiscreteUniform) = 0.0

function kurtosis(d::DiscreteUniform)
    n2 = span(d)^2
    -1.2 * (n2 + 1.0) / (n2 - 1.0)
end

entropy(d::DiscreteUniform) = log(span(d))

mode(d::DiscreteUniform) = d.a
modes(d::DiscreteUniform) = [d.a:d.b]

### Evaluation

pdf(d::DiscreteUniform, x::Real) = insupport(d, x) ? d.pv : zero(d.pv)
logpdf(d::DiscreteUniform, x::Real) = log(pdf(d, x))

function cdf(d::DiscreteUniform, x::Int)
    a = d.a
    result = (x - a + 1) * d.pv
    return if x < a
        zero(result)
    elseif x >= d.b
        one(result)
    else
        result
    end
end

quantile(d::DiscreteUniform, p::Real) = iszero(p) ? d.a : d.a - 1 + ceil(Int, p * span(d))

function mgf(d::DiscreteUniform, t::Real)
    a, b = d.a, d.b
    u = b - a + 1
    result = (exp(t*a) * expm1(t*u)) / (u*expm1(t))
    return iszero(t) ? one(result) : result
end

function cf(d::DiscreteUniform, t::Real)
    a, b = d.a, d.b
    u = b - a + 1
    result = (im*cos(t*(a+b)/2) + sin(t*(a-b-1)/2)) / (u*sin(t/2))
    return iszero(t) ? one(result) : result
end

### Sampling

rand(rng::AbstractRNG, d::DiscreteUniform) = rand(rng, d.a:d.b)

# Fit model

function fit_mle(::Type{DiscreteUniform}, x::AbstractArray{<:Real})
    if isempty(x)
        throw(ArgumentError("data set must be non-empty."))
    end
    return DiscreteUniform(extrema(x)...)
end


text:
function maximum_weight_maximal_matching_hungarian(g::Graph,
          w::AbstractMatrix{T}=default_weights(g)) where {T <: Real}
  edge_list = collect(edges(g))
  n = nv(g)

  # Determine the bipartition of the graph.
  bipartition = bipartite_map(g)
  if length(bipartition) != n # Equivalent to !is_bipartite(g), but reuses the results of the previous function call.
    error("The Hungarian algorithm only works for bipartite graphs; otherwise, prefer the Blossom algorithm (not yet available in GraphsMatching")
  end
  n_first = count(bipartition .== 1)
  n_second = count(bipartition .== 2)

  to_bipartition_1 = [count(bipartition[1:i] .== 1) for i in 1:n]
  to_bipartition_2 = [count(bipartition[1:i] .== 2) for i in 1:n]

  # hungarian() minimises the total cost, while this function is supposed to maximise the total weights.
  wDual = maximum(w) .- w

  # Remove weights that are not in the graph (Hungarian.jl considers all weights that are not missing values as real edges).
  # Assume w is symmetric, so that the weight of matching i->j is the same as the one for j->i.
  weights = Matrix{Union{Missing, T}}(missing, n_first, n_second)

  for i in 1:n
    for j in 1:n
      if Edge(i, j) ∈ edge_list || Edge(j, i) ∈ edge_list
        if bipartition[i] == 1 # and bipartition[j] == 2
          idx_first = to_bipartition_1[i]
          idx_second = to_bipartition_2[j]
        else # bipartition[i] == 2 and bipartition[j] == 1
          idx_first = to_bipartition_1[j]
          idx_second = to_bipartition_2[i]
        end

        weight_to_add = (Edge(i, j) ∈ edge_list) ? wDual[i, j] : wDual[j, i]

        weights[idx_first, idx_second] = weight_to_add
      end
    end
  end

  # Run the Hungarian algorithm.
  assignment, _ = hungarian(weights)

  # Convert the output format to match LGMatching's.
  pairs = Tuple{Int, Int}[]
  mate = fill(-1, n) # Initialise to unmatched.
  for i in eachindex(assignment)
    if assignment[i] != 0 # If matched:
      original_i = findfirst(to_bipartition_1 .== i)
      original_j = findfirst(to_bipartition_2 .== assignment[i])

      mate[original_i] = original_j
      mate[original_j] = original_i

      push!(pairs, (original_i, original_j))
    end
  end

  # Compute the cost for this matching (as weights had to be changed for Hungarian.jl, the one returned by hungarian() makes no sense).
  cost = sum(w[p[1], p[2]] for p in pairs)

  # Return the result.
  return MatchingResult(cost, mate)
end


text:
using Dates
using ArgParse
using Printf
using Random
using GaussianProcesses

# before_all = Dates.now()

include("../../btg.jl")

###
###  Manually enter parsed args entries so we can run this file as a script
###
##
parsed_args = Dict()
push!(parsed_args, "randseed" => 1234)
push!(parsed_args, "lmax" => 0.5)
push!(parsed_args, "lmin" => 0.2)
push!(parsed_args, "noiselevel" => 1)
push!(parsed_args, "ntest" => 200)
push!(parsed_args, "ntrain" => 80)
push!(parsed_args, "p" => 1)

randseed = parsed_args["randseed"]; rng = MersenneTwister(randseed)
Random.seed!(randseed);

function covariate_fun(x, p)
    n = size(x, 1)
    d = size(x, 2)
    if p == 1
        return ones(n, 1)
    elseif p == 1 + d
        return hcat(ones(n), x)
    else
        throw(ArgumentError("Only support constant or linear convariate."))
    end
end

true_fun(x, s) = sin(x) + s*randn() + 3
true_fun_noise(x) = max(true_fun(x, 1e-1), 0)^1/3
true_fun_0(x) = true_fun(x, 0.)^1/3

# training set
n_train = parsed_args["ntrain"]
data = range(-pi, stop=pi, length=n_train)
#data = sort((rand(Float64, (n_train, )) .- 0.5)*2*pi)
target = (sin.(data) .+ 0.1 * parsed_args["noiselevel"].* randn(rng, n_train) .+ 3).^(1/3)
x = reshape(data, n_train, 1)
Fx = covariate_fun(x, parsed_args["p"])
max_train = maximum(target)
y = target ./ max_train
trainingData0 = trainingData(x, Fx, y)
d = getDimension(trainingData0); n = getNumPts(trainingData0); p = getCovDimension(trainingData0)
# testing set
n_test = parsed_args["ntest"]
x_test = range(-pi, stop=pi, length=n_test)
#x_test = (rand(Float64, (n_test, )) .- 0.5)*2*pi
x_test = sort(reshape(x_test, n_test, 1), dims=1)
Fx_test = covariate_fun(x_test, parsed_args["p"])
y_test_true = (sin.(x_test) .+ 3).^(1/3)
elapsedmin = 0

#parameter setting
# myquadtype = parsed_args["sparse"] ? ["SparseCarlo", "SparseCarlo"] : ["QuasiMonteCarlo", "QuasiMonteCarlo"]
myquadtype = ["Gaussian", "Gaussian"]
rangeλ = [-20.0 20.0]
#rangeλ = reshape([0.1], 1, 1)
lmin = parsed_args["lmin"]
lmax = parsed_args["lmax"]
#rangeθ = [1/lmax^2 1/lmin^2]
rangeθ = reshape([1.0], 1, 1)
@info "rangeθ, rangel:" rangeθ, [lmin lmax]
# rangeθ = [0.111 25]
# build btg model

function load_synthetic_btg()
    return btg(trainingData0, rangeθ, rangeλ; priorθ = inverseUniform(rangeθ), quadtype = myquadtype)
end

function getseed()
    return randseed
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
import SGP4
using Base.Test

line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  4753"
line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.82419157413667"

wgs72 = SGP4.GravityModel("wgs72")
satellite = SGP4.twoline2rv(line1, line2, wgs72)

@test satellite[:satnum] == 5

t = Dates.DateTime(2000, 6, 29, 12, 50, 19)
(position, velocity) = SGP4.propagate( satellite, 2000, 6, 29, 12, 50, 19)
(position2, velocity2) = SGP4.propagate( satellite, t )

@test satellite[:error] == 0

@test_approx_eq_eps position[1] position2[1] eps()
@test_approx_eq_eps position[2] position2[2] eps()
@test_approx_eq_eps position[3] position2[3] eps()

@test_approx_eq_eps velocity[1] velocity2[1] eps()
@test_approx_eq_eps velocity[2] velocity2[2] eps()
@test_approx_eq_eps velocity[3] velocity2[3] eps()

@test_approx_eq_eps position[1] 5576.056952 1e-6
@test_approx_eq_eps position[2] -3999.371134 1e-6
@test_approx_eq_eps position[3] -1521.957159 1e-6

@test_approx_eq_eps velocity[1] 4.772627 1e-6
@test_approx_eq_eps velocity[2] 5.119817 1e-6
@test_approx_eq_eps velocity[3] 4.275553 1e-6

sats = [satellite; satellite; satellite]

(positions, velocities) = SGP4.propagate( sats, t )

@test_approx_eq_eps positions[1,1] positions[1,3] eps()
@test_approx_eq_eps positions[2,1] positions[2,3] eps()
@test_approx_eq_eps positions[3,1] positions[3,3] eps()

@test_approx_eq_eps velocities[1,1] velocities[1,3] eps()
@test_approx_eq_eps velocities[2,1] velocities[2,3] eps()
@test_approx_eq_eps velocities[3,1] velocities[3,3] eps()


text:
"""
    classnames() -> Vector{Int}

Return the 10 digits for the SVHN classes as a vector of integers.
"""
classnames() = CLASSES

for (FUN, PATH, COUNT, DESC) in (
        (:traintensor, TRAINDATA, 73257,  "training"),
        (:testtensor,  TESTDATA,  26032,  "test"),
        (:extratensor, EXTRADATA, 531131, "extra training"),
       )
    FUN_STR = string(FUN)
    @eval begin
        """
            $($FUN_STR)([T = N0f8], [indices]; [dir]) -> Array{T}

        Return the SVHN **$($DESC)** images corresponding to the
        given `indices` as a multi-dimensional array of eltype
        `T`.

        The image(s) is/are returned in the native vertical-major
        memory layout as a single numeric array. If `T <:
        Integer`, then all values will be within `0` and `255`,
        otherwise the values are scaled to be between `0` and
        `1`.

        If the parameter `indices` is omitted or an
        `AbstractVector`, the images are returned as a 4D array
        (i.e. a `Array{T,4}`), in which the first dimension
        corresponds to the pixel *columns* (y) of the image, the
        second dimension to the pixel *rows* (x) of the image,
        the third dimension the RGB color channels, and the
        fourth dimension denotes the index of the image.

        ```julia-repl
        julia> SVHN2.$($FUN_STR)() # load all $($DESC) images
        32×32×3×$($COUNT) Array{N0f8,4}:
        [...]

        julia> SVHN.$($FUN_STR)(Float32, 1:3) # first three images as Float32
        32×32×3×3 Array{Float32,4}:
        [...]
        ```

        If `indices` is an `Integer`, the single image is
        returned as `Array{T,3}` in vertical-major layout, which
        means that the first dimension denotes the pixel
        *columns* (y), the second dimension denotes the pixel
        *rows* (x), and the third dimension the RGB color
        channels of the image.

        ```julia-repl
        julia> SVHN2.$($FUN_STR)(1) # load first $($DESC) image
        32×32×3 Array{N0f8,3}:
        [...]
        ```

        As mentioned above, the color channel is encoded in the
        third dimension. You can use the utility function
        [`convert2image`](@ref) to convert an SVHN array into a
        Julia image with the appropriate `RGB` eltype.

        ```julia-repl
        julia> SVHN2.convert2image(SVHN2.$($FUN_STR)(1))
        32×32 Array{RGB{N0f8},2}:
        [...]
        ```

        $(download_docstring("SVHN2", DEPNAME))
        """
        function ($FUN)(args...; dir = nothing)
            ($FUN)(N0f8, args...; dir = dir)
        end

        function ($FUN)(::Type{T}; dir = nothing) where T
            path = datafile(DEPNAME, $PATH, dir)
            images = matopen(io->read(io, "X"), path)::Array{UInt8,4}
            images = permutedims(images, (2, 1, 3, 4))
            bytes_to_type(T, images)
        end

        function ($FUN)(::Type{T}, indices; dir = nothing) where T
            images = ($FUN)(T, dir = dir)#::Array{T,4}
            images[:,:,:,indices]
        end
    end
end

for (FUN, PATH, COUNT, DESC) in (
        (:trainlabels, TRAINDATA, 73257,  "training"),
        (:testlabels,  TESTDATA,  26032,  "test"),
        (:extralabels, EXTRADATA, 531131, "extra training"),
       )
    FUN_STR = string(FUN)
    @eval begin
        """
            $($FUN_STR)([indices]; [dir])

        Returns the SVHN **$($DESC)** labels corresponding to
        the given `indices` as an `Int` or `Vector{Int}`. The
        values of the labels denote the zero-based class-index
        that they represent (see [`SVHN2.classnames`](@ref) for
        the corresponding names). If `indices` is omitted, all
        labels are returned.

        ```julia-repl
        julia> SVHN2.$($FUN_STR)() # full $($DESC) set
        $($COUNT)-element Array{Int64,1}:
        [...]

        julia> SVHN2.$($FUN_STR)(1:3) # first three labels
        3-element Array{Int64,1}:
        [...]

        julia> SVHN2.$($FUN_STR)(1) # first label
        [...]

        julia> SVHN2.classnames()[SVHN2.$($FUN_STR)(1)] # corresponding class
        [...]
        ```

        $(download_docstring("SVHN2", DEPNAME))
        """
        function ($FUN)(; dir = nothing)
            path = datafile(DEPNAME, $PATH, dir)
            labels = matopen(io->read(io, "y"), path)
            Vector{Int}(vec(labels))::Vector{Int}
        end

        function ($FUN)(indices; dir = nothing)
            labels = ($FUN)(dir = dir)::Vector{Int}
            labels[indices]
        end
    end
end

for (FUN, PATH, DESC) in (
        (:traindata, TRAINDATA, "trainset"),
        (:testdata,  TESTDATA,  "testset"),
        (:extradata, EXTRADATA, "extra trainset"),
       )
    FUN_STR = string(FUN)
    @eval begin
        """
            $($FUN_STR)([T = N0f8], [indices]; [dir]) -> images, labels

        Returns the SVHN **$($DESC)** corresponding to the given
        `indices` as a two-element tuple. If `indices` is omitted
        the full $($DESC) is returned. The first element of the
        return values will be the images as a multi-dimensional
        array, and the second element the corresponding labels as
        integers.

        The image(s) is/are returned in the native vertical-major
        memory layout as a single numeric array of eltype `T`. If
        `T <: Integer`, then all values will be within `0` and
        `255`, otherwise the values are scaled to be between `0`
        and `1`. You can use the utility function
        [`convert2image`](@ref) to convert an SVHN array into a
        Julia image with the appropriate `RGB` eltype. The
        integer values of the labels correspond 1-to-1 the digit
        that they represent with the exception of 0 which is
        encoded as `10`.

        Note that because of the nature of how the dataset is
        stored on disk, `SVHN2.$($FUN_STR)` will always load the
        full $($DESC), regardless of which observations are
        requested. In the case `indices` are provided by the
        user, it will simply result in a sub-setting. This option
        is just provided for convenience.

        ```julia
        images, labels = SVHN2.$($FUN_STR)() # full dataset
        images, labels = SVHN2.$($FUN_STR)(2) # only second observation
        images, labels = SVHN2.$($FUN_STR)(dir="./SVHN") # custom folder
        ```

        $(download_docstring("SVHN2", DEPNAME))
        """
        function ($FUN)(args...; dir = nothing)
            ($FUN)(N0f8, args...; dir = dir)
        end

        function ($FUN)(::Type{T}; dir = nothing) where T
            path = datafile(DEPNAME, $PATH, dir)
            vars = matread(path)
            images = vars["X"]::Array{UInt8,4}
            labels = vars["y"]
            images = permutedims(images, (2, 1, 3, 4))
            bytes_to_type(T, images), Vector{Int}(vec(labels))::Vector{Int}
        end

        function ($FUN)(::Type{T}, indices; dir = nothing) where T
            images, labels = ($FUN)(T, dir = dir)
            images[:,:,:,indices], labels[indices]
        end
    end
end


text:
using Plots: string
using Base: Float16, Float64
include("PolygonBase.jl")
include("ESBase.jl")
include("silhouette&pieces.jl")

levy(X::Array) = sin(3π*X[1])^2 + (X[1] - 1)^2 * (1 + sin(3π*X[2])^2) + (X[2] - 1)^2 * (1 + sin(2π * X[2])^2)
levy(x, y) = levy([x, y])

rastrigin(x, y) = 10 * 2 + x^2 - 10cos(2π*x) + y^2 - 10cos(2π*y)
rastrigin(X::Array) = rastrigin(X[1], X[2])

eggholder(x, y) = -(100*y + 47)*sin(sqrt(abs(50x + 100y + 47))) - 100x*sin(sqrt(abs(100x - 100y - 47)))
eggholder(X::Array) = eggholder(X[1], X[2])

# X[1つめのピースのx座標, 1つめのピースのy座標, １つめの回転角Θ, 2つめのx座標, ...]
function loss_poly(X::Array, loss_args::Array{Float64, 1})
    global pices, silhouette

    # silhouetteとpieceの共通部分の面積を求める
    tmp_p = move(pieces[1], X[1], X[2])
    rotate!(tmp_p, X[3] * 2π)
    unioned = MYPolygon2LibGEOS(tmp_p)
    sum_pieces = LibGEOS.area(MYPolygon2LibGEOS(tmp_p))

    for i in 2:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        unioned = LibGEOS.union(unioned, MYPolygon2LibGEOS(tmp_p))
        sum_pieces += MYPolygon2LibGEOS(tmp_p) |> LibGEOS.area
    end
    tmp = LibGEOS.intersection(unioned, MYPolygon2LibGEOS(silhouette))
    A = LibGEOS.area(tmp) / LibGEOS.area(MYPolygon2LibGEOS(silhouette))

    # 各2pieceの共通部分を求める
    inter = 1.0 - LibGEOS.area(unioned) / sum_pieces
    outer = 1.0 - LibGEOS.area(tmp)/sum_pieces

    # 頂点の誤差成分 (シルエットベース 緩い)
    Δv = 0.0
    for v_s in 1:silhouette.n
        Δmin = 100.0
        for j in 1:length(pieces)
            tmp_p = move(pieces[j], X[3j - 2], X[3j - 1])
            rotate!(tmp_p, X[3j] * 2π)
            for v_p in 1:tmp_p.n
                Δmin = min(Δmin, sum((silhouette.vertexes[:,v_s] - tmp_p.vertexes[:, v_p]).^2))
            end
        end
        Δv += Δmin / silhouette.n
        # Δv = max(Δv, Δmin)
    end

    # 頂点の誤差成分 (ピースベース 厳しい)
    """
    Δv = 0.0
    for i in 1:length(pieces)
        p = move(pieces[i], X[3i - 2], X[3i - 1])
        rotate!(p, X[3i] * 2π)
        Δmax = 0.0
        # 頂点ごとに誤差を見る
        for v_p in 1:p.n
            Δmin = 100.0
            # シルエットについて見る
            for v_s in 1:silhouette.n
                Δmin = min(Δmin, sum((silhouette.vertexes[:,v_s] - p.vertexes[:, v_p]).^2))
            end
            # その他のピースについて見る
            for j in 1:length(pieces)
                if j == i continue end
                q = move(pieces[j], X[3j - 2], X[3j - 1])
                rotate!(p, X[3j] * 2π)
                for v_q in 1:q.n
                    Δmin = min(Δmin, sum((q.vertexes[:,v_q] - p.vertexes[:, v_p]).^2))
                end
            end
            Δmax = max(Δmax, Δmin)
        end
        # Δv += 2tanh(2Δmin) / (π*silhouette.n)
        # Δv = max(Δv, Δmin)
        Δv += Δmax
    end
    """

    # 辺の誤差成分
    E = 0.0
    for i in 1:length(pieces)
        P = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(P, X[3i] * 2π)

        cosθ = 0.0

        # シルエットの各頂点と距離が最小な点を探す
        min_d = 100.0
        for p in 0:P.n-1, s in 0:silhouette.n-1
            d = sum((P.vertexes[:, p + 1] - silhouette.vertexes[:, s + 1]).^2)
            # 最小の大きさの角度を探す
            if d < min_d
                eP = P.vertexes[:, mod(p, P.n) + 1] - P.vertexes[:, mod(p-1, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s, silhouette.n) + 1] - silhouette.vertexes[:, mod(s-1, silhouette.n) + 1]
                cosθ1 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p, P.n) + 1] - P.vertexes[:, mod(p-1, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s+1, silhouette.n) + 1] - silhouette.vertexes[:, mod(s, silhouette.n) + 1]
                cosθ2 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p+1, P.n) + 1] - P.vertexes[:, mod(p, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s, silhouette.n) + 1] - silhouette.vertexes[:, mod(s-1, silhouette.n) + 1]
                cosθ3 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p+1, P.n) + 1] - P.vertexes[:, mod(p, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s+1, silhouette.n) + 1] - silhouette.vertexes[:, mod(s, silhouette.n) + 1]
                cosθ4 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                cosθ = max(cosθ1, cosθ2, cosθ3, cosθ4)
            end
        end

        E += √(1 - cosθ^2) / length(pieces)
    end

    return -loss_args[1] * A + loss_args[2] * Δv + loss_args[3] * outer + loss_args[4] * E
end

silhouette = hexagon_5
pieces = [tri_m, tri_s, tri_s, square_s, parallelogram]

# 初期化
cmaes = init_CMAES(zeros(3 * length(pieces)), 1.0, 0)
seed = MersenneTwister()
max_gen = 168

println(seed)

fitnesses_ave = Array{Float64, 1}()
fitnesses_max = Array{Float64, 1}()

loss_args = [100.0, 20.0, 50.0, 20.0]

anim = @animate for gen in 1:max_gen
    global fitnesses_ave, fitnesses_max, seed, cmaes, loss_args
    local X, fitnesses

    # 個体生成
    X = samplePopulation(cmaes, rng=seed)

    # 回転角のパラメータを -π ~ π までに正規化
    for j in 1:size(X)[2]
        for i in 3:3:size(X)[1]
            X[i, j] = mod(X[i, j], 1.0)
        end
    end

    # 評価値
    loss(x) = loss_poly(x, loss_args)
    fitnesses = get_fitness(X, loss)

    # 表示用
    best_arg = argmin(fitnesses)
    best_pieces = Array{MYPolygon, 1}()
    for i in 1:length(pieces)
        todisplay = move(pieces[i], X[3i-2, best_arg], X[3i-1, best_arg])
        rotate!(todisplay, X[3i, best_arg] * 2π)
        push!(best_pieces, todisplay)
    end
    display(silhouette, best_pieces...)

    # 評価用
    push!(fitnesses_ave, sum(fitnesses)/size(X)[2])
    push!(fitnesses_max, minimum(fitnesses))

    # 更新
    update!(cmaes, X, fitnesses, gen % (max_gen ÷ 2))

    # 再出発
    if false # 2gen == max_gen
        # println("Searching in local...")
        loss_args = [100.0, 0.0, 100.0, 0.0]
        best_arg = argmin(fitnesses)
        cmaes = init_CMAES(X[:, best_arg], 1.0, 128)
    end
end

gif(anim, "best_pieces.gif", fps=10)

pl = plot(1:max_gen, -fitnesses_max, lab="Average", xaxis="generation", yaxis="fitness [%]", legend = false, w=2.5)
# plot!(pl, 1:max_gen, -fitnesses_max, line=:dash, lab="Maximum", ylim=(minimum(-fitnesses_ave), 110), yticks = 100:-20:minimum(-fitnesses_ave))
plot!(pl, xtickfontsize=13, ytickfontsize=13, xguidefontsize=13, yguidefontsize=13)

savefig(pl, "fitnesses.png")
Base.display(pl)

max_stars_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_issues_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_forks_repo_name: toshiki-kosen/SolvingTangramPuzzle


text:
const testmodel = Catwalk.DefaultDispatchCostModel(
    skip                = 1,
    static_dispatch     = 10,
    dynamic_dispatch    = 100,
)

const testfreqs = Catwalk.DataTypeFrequencies()
for i =1:3 Catwalk.increment!(testfreqs, Int) end
for i =1:111 Catwalk.increment!(testfreqs, Float16) end
for i =1:5 Catwalk.increment!(testfreqs, Float32) end
for i =1:7 Catwalk.increment!(testfreqs, Float64) end

@testset "Cost Model" begin
    @test Catwalk.costof(values(testfreqs), Catwalk.encode(Float16, Int, Float32), testmodel) ==
    (
        3   * (2 * testmodel.skip + testmodel.static_dispatch) +
        111 * (1 * testmodel.skip + testmodel.static_dispatch) +
        5   * (3 * testmodel.skip + testmodel.static_dispatch) +
        7   * (3 + testmodel.dynamic_dispatch)
    )
end


text:
"""
    common.jl

Description:
    Includes all of the unsupervised ARTMAP modules common code.
"""

# -------------------------------------------
# Methods
# -------------------------------------------

"""
    train!(art::ARTMAP, x::RealMatrix, y::IntegerVector, preprocessed::Bool=false)

Train the ARTMAP model on a batch of data 'x' with supervisory labels 'y.'

# Arguments
- `art::ARTMAP`: the supervised ARTMAP model to train.
- `x::RealMatrix`: the 2-D dataset containing columns of samples with rows of features.
- `y::IntegerVector`: labels for supervisory training.
- `preprocessed::Bool=false`: flag, if the data has already been complement coded or not.
"""
function train!(art::ARTMAP, x::RealMatrix, y::IntegerVector, preprocessed::Bool=false)
    # Show a message if display is on
    art.opts.display && @info "Training $(typeof(art))"

    # Data information and setup
    n_samples = length(y)

    # Run the batch initialization procedure
    x = init_train!(x, art, preprocessed)

    # Initialize the output vector
    y_hat = zeros(Int, n_samples)
    # Learn until the stopping conditions
    art.epoch = 0
    while true
        # Increment the epoch and get the iterator
        art.epoch += 1
        iter = get_iterator(art.opts, x)
        for i = iter
            # Update the iterator if necessary
            update_iter(art, iter, i)
            # Grab the sample slice
            # sample = get_sample(x, i)
            sample = x[:, i]
            label = y[i]
            # Train upon the sample and label
            y_hat[i] = train!(art, sample, label, preprocessed=true)
        end

        # Check stopping conditions
        if stopping_conditions(art)
            break
        end
    end
    return y_hat
end # train!(art::ARTMAP, x::RealMatrix, y::IntegerVector, preprocessed::Bool=false)

# -------------------------------------------
# Common Documentation
# -------------------------------------------

@doc raw"""
    train!(art::ARTMAP, x::RealVector, y::Integer ; preprocessed::Bool=false)

Train the supervised ARTMAP model on a single sample of features 'x' with supervisory label 'y'.

# Arguments
- `art::ARTMAP`: the supervised ART model to train.
- `x::RealVector`: the single sample feature vector to train upon.
- `y::Integer`: the label for supervisory training.
- `preprocessed::Bool=false`: optional, flag if the data has already been complement coded or not.
"""
train!(art::ARTMAP, x::RealVector, y::Integer ; preprocessed::Bool=false)


text:
@testset "CURAND" begin

using CuArrays.CURAND

CURAND.seed!()
CURAND.seed!(1)
CURAND.seed!(1, 0)

rng = CURAND.generator()
CURAND.seed!(rng)
CURAND.seed!(rng, 1)
CURAND.seed!(rng, 1, 0)

# NOTE: tests should cover both pow2 and non-pow2 dims

# in-place
for (f,T) in ((rand!,Float32),
              (randn!,Float32),
              (rand_logn!,Float32),
              (rand_poisson!,Cuint)),
    d in (2, (2,2), (2,2,2), 3, (3,3), (3,3,3))
    A = CuArray{T}(undef, d)
    f(A)
end

# out-of-place, with implicit type
for (f,T) in ((CuArrays.rand,Float32), (CuArrays.randn,Float32),
              (CuArrays.rand_logn,Float32), (CuArrays.rand_poisson,Cuint),
              (rand,Float64), (randn,Float64)),
    args in ((2,), (2, 2), (3,), (3, 3))
    A = f(args...)
    @test eltype(A) == T
end

# out-of-place, with type specified
for (f,T) in ((CuArrays.rand,Float32), (CuArrays.randn,Float32), (CuArrays.rand_logn,Float32),
              (CuArrays.rand,Float64), (CuArrays.randn,Float64), (CuArrays.rand_logn,Float64),
              (CuArrays.rand_poisson,Cuint),
              (rand,Float32), (randn,Float32),
              (rand,Float64), (randn,Float64)),
    args in ((T, 2), (T, 2, 2), (T, (2, 2)), (T, 3), (T, 3, 3), (T, (3, 3)))
    A = f(args...)
    @test eltype(A) == T
end

# unsupported types that fall back to GPUArrays
for (f,T) in ((CuArrays.rand,Int64),),
    args in ((T, 2), (T, 2, 2), (T, (2, 2)), (T, 3), (T, 3, 3), (T, (3, 3)))
    A = f(args...)
    @test eltype(A) == T
end
for (f,T) in ((rand!,Int64),),
    d in (2, (2,2), (2,2,2), 3, (3,3), (3,3,3))
    A = CuArray{T}(undef, d)
    f(A)
end

@test_throws ErrorException randn!(CuArray{Cuint}(undef, 10))
@test_throws ErrorException rand_logn!(CuArray{Cuint}(undef, 10))
@test_throws ErrorException rand_poisson!(CuArray{Float64}(undef, 10))

end


text:
function DayAheadModel(data::DayAheadData, scenarios::Vector{<:DayAheadScenario}; kw...)
    return StochasticHydroModel("Day-Ahead", HydroModels.Day(), data, DayAheadModelDef, scenarios; kw...)
end

function DayAheadModel(data::DayAheadData, sampler::RecurrentDayAheadSampler, n::Integer; kw...)
    return StochasticHydroModel("Day-Ahead", HydroModels.Day(), data, DayAheadModelDef, sampler, n; kw...)
end

function DayAheadModelDef(horizon::Horizon, data::DayAheadData, indices::DayAheadIndices)
    stochasticmodel = @stochastic_model begin
        # First stage
        # ========================================================
        @stage 1 begin
            @parameters begin
                horizon = horizon
                indices = indices
                data = data
            end
            @unpack hours, plants, bids, blockbids, blocks, hours_per_block = indices
            @unpack hydrodata, regulations, use_blockbids = data
            active_blocks(t) = findall(A->in(t,A), hours_per_block)
            # Variables
            # ========================================================
            @decision(model, xᴵ[t in hours] >= 0)
            @decision(model, xᴰ[i in bids, t = hours] >= 0)
            if use_blockbids
                @decision(model, 0 <= xᴮ[i = blockbids, b = blocks] <= regulations.blocklimit)
            end
            # Objectives
            # ========================================================
            # Dummy first-stage objective in order to use Max sense
            @objective(model, Max, 0)
            # Constraints
            # ========================================================
            # Increasing bid curve
            @constraint(model, bidcurve[i in bids[1:end-1], t in hours],
                xᴰ[i,t] <= xᴰ[i+1,t]
            )
            # Maximal bids
            if use_blockbids
                @constraint(model, maxhourlybids[t in hours],
                    xᴵ[t] + xᴰ[bids[end],t] + sum(xᴮ[i,b]
                        for i in blockbids for b in active_blocks(t))
                            <= 2*sum(hydrodata[p].H̄ for p in plants)
                )
            else
                @constraint(model, maxhourlybids[t in hours],
                    xᴵ[t] + xᴰ[bids[end],t] <= 2*sum(hydrodata[p].H̄ for p in plants)
                )
            end
        end
        # Second stage
        # =======================================================
        @stage 2 begin
            @parameters begin
                horizon = horizon
                indices = indices
                data = data
            end
            @unpack hours, plants, segments, blocks, hours_per_block = indices
            @unpack hydrodata, water_value, regulations, intraday_trading, penalty_percentage, use_blockbids, bidlevels = data
            α = penalty_percentage
            β = α == 0.0 ? 0.0 : 1/α
            @uncertain ρ Q̃ from ξ::DayAheadScenario
            V = local_inflows(Q̃, hydrodata.Qu)
            # Auxilliary functions
            ih(t) = begin
                idx = findlast(bidlevels[t] .<= ρ[t])
                return idx == nothing ? -1 : idx
            end
            function blockbidlevel(b)
                return mean([bidlevels[i][2:end-1] for i in hours_per_block[b]])
            end
            ib(b) = begin
                idx = findlast(blockbidlevel(b) .<= mean(ρ[hours_per_block[b]]))
                return idx == nothing ? -1 : idx
            end
            interpolate(ρ, bidlevels, xᴰ, t) = begin
                lower = ((ρ[t] - bidlevels[t][ih(t)])/(bidlevels[t][ih(t)+1]-bidlevels[t][ih(t)]))*xᴰ[ih(t)+1,t]
                upper = ((bidlevels[t][ih(t)+1]-ρ[t])/(bidlevels[t][ih(t)+1]-bidlevels[t][ih(t)]))*xᴰ[ih(t),t]
                return lower + upper
            end
            active_blocks(t) = findall(A->in(t,A), hours_per_block)
            # Variables
            # =======================================================
            # -------------------------------------------------------
            @variable(model, yᴴ[t in hours] >= 0)
            if use_blockbids
                @variable(model, yᴮ[b in blocks] >= 0)
            end
            if intraday_trading
                @variable(model, y⁺[t in hours] >= 0)
                @variable(model, y⁻[t in hours] >= 0)
            end
            @variable(model, 0 <= Q[p in plants, s in segments, t in hours] <= Q̄(hydrodata, p, s))
            @variable(model, S[p in plants, t in hours] >= 0)
            @variable(model, 0 <= M[p in plants, t in hours] <= hydrodata[p].M̄)
            @variable(model, W[i in 1:nindices(water_value)])
            @variable(model, H[t in hours] >= 0)
            @variable(model, Qf[p in plants, t in hours] >= 0)
            @variable(model, Sf[p in plants, t in hours] >= 0)

            # Objectives
            # ========================================================
            # Net profit
            if use_blockbids
                @expression(model, net_profit,
                    sum((ρ[t]-regulations.dayaheadfee)*yᴴ[t]
                        for t in hours)
                    + sum(length(hours_per_block[b])*(mean(ρ[hours_per_block[b]])-regulations.dayaheadfee)*yᴮ[b]
                        for b in blocks))
            else
                @expression(model, net_profit,
                    sum((ρ[t]-regulations.dayaheadfee)*yᴴ[t]
                        for t in hours))
            end
            # Intraday
            if intraday_trading
                @expression(model, intraday,
                            sum((penalty(ξ,t) + regulations.intradayfee)*y⁺[t] - (reward(ξ,t) - regulations.intradayfee)*y⁻[t]
                                for t in hours))
            end
            # Value of stored water
            @expression(model, value_of_stored_water, -sum(W[i] for i in 1:nindices(water_value)))
            # Define objective
            if intraday_trading
                @objective(model, Max, net_profit - intraday + value_of_stored_water)
            else
                @objective(model, Max, net_profit + value_of_stored_water)
            end

            # Constraints
            # ========================================================
            # Bid-dispatch links
            @constraint(model, hourlybids[t in hours],
                yᴴ[t] == interpolate(ρ, bidlevels, xᴰ, t) + xᴵ[t]
            )
            if use_blockbids
                @constraint(model, bidblocks[b in blocks],
                    yᴮ[b] == sum(xᴮ[j,b]
                        for j in 1:ib(b)))
            end

            # Hydrological balance
            @constraint(model, hydro_constraints[p in plants, t in hours],
                # Previous reservoir content
                M[p,t] == (t > 1 ? M[p,t-1] : hydrodata[p].M₀)
                # Inflow
                + sum(Qf[i,t] for i in intersect(hydrodata.Qu[p],plants))
                + sum(Sf[i,t] for i in intersect(hydrodata.Su[p],plants))
                # Local inflow
                + V[p]
                # Outflow
                - sum(Q[p,s,t] for s in segments)
                - S[p,t]
            )
            # Production
            @constraint(model, production[t in hours],
                        H[t] == sum(marginal_production(Resolution(1), μ(hydrodata, p, s)) * Q[p,s,t]
                                    for p in plants, s in segments)
                        )
            # Load balance
            if intraday_trading
                if use_blockbids
                    @constraint(model, loadbalance[t in hours],
                        yᴴ[t] + sum(yᴮ[b] for b in active_blocks(t)) - H[t] == y⁺[t] - y⁻[t]
                    )
                else
                    @constraint(model, loadbalance[t in hours],
                        yᴴ[t] - H[t] == y⁺[t] - y⁻[t]
                    )
                end
            else
                if use_blockbids
                    @constraint(model, loadbalance[t in hours],
                        yᴴ[t] + sum(yᴮ[b] for b in active_blocks(t)) == H[t]
                    )
                else
                    @constraint(model, loadbalance[t in hours],
                        yᴴ[t] == H[t]
                    )
                end
            end
            # Water flow: Discharge + Spillage
            @constraint(model, discharge_flow_time[p in plants, t in hours],
                        Qf[p,t] == (t - water_flow_time(Resolution(1), hydrodata[p].Rq) > 0 ?
                        overflow(Resolution(1), hydrodata[p].Rq) * sum(Q[p,s,t-water_flow_time(Resolution(1), hydrodata[p].Rq)]
                                                                    for s in segments) : 0.0)
                        + (t - water_flow_time(Resolution(1), hydrodata[p].Rq) > 1 ?
                        historic_flow(Resolution(1), hydrodata[p].Rq)*sum(Q[p,s,t-(water_flow_time(Resolution(1), hydrodata[p].Rq)+1)]
                                                                       for s in segments) : 0.0))
            @constraint(model, spillage_flow_time[p in plants, t in hours],
                        Sf[p,t] == (t - water_flow_time(Resolution(1), hydrodata[p].Rs) > 0 ?
                        overflow(Resolution(1), hydrodata[p].Rs) * S[p,t-water_flow_time(Resolution(1), hydrodata[p].Rs)] : 0.0)
                        + (t - water_flow_time(Resolution(1), hydrodata[p].Rs) > 1 ?
                        historic_flow(Resolution(1), hydrodata[p].Rs)*S[p,t-(water_flow_time(Resolution(1), hydrodata[p].Rq)+1)] : 0.0))
            # Water value
            @constraint(model, water_value_approximation[c in 1:ncuts(water_value)],
                        sum(water_value[c][p]*M[p,num_hours(horizon)]
                            for p in plants)
                        + sum(W[i]
                              for i in cut_indices(water_value[c])) >= cut_lb(water_value[c]))
        end
    end
    return stochasticmodel
end


text:
function execute_code(code)
    acc = 0
    curr = 1
    seen = Set{Int}()
    while curr != length(code) + 1
        if curr in seen
            return nothing
        end
        push!(seen, curr)
        op, arg = code[curr]
        if op == "acc"
            curr += 1
            acc += arg
        elseif op == "jmp"
            curr += arg
        else
            curr += 1
        end
    end
    return acc
end

function run(s)
    # Your code here
    code = readlines(IOBuffer(s))
    code = [(line[1:3], parse(Int, line[5:end])) for line in code]
    for (idx, (op, arg)) in enumerate(code)
        if op == "acc"
            continue
        elseif op == "jmp"
            code[idx] = ("nop", arg)
        else
            code[idx] = ("jmp", arg)
        end
        acc = execute_code(code)
        if acc === nothing
            code[idx] = (op, arg)
        else
            return acc
        end
    end
    return 0
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
using SafeTestsets

@safetestset "filtered chain" begin
    include("filteredchain.jl")
end
@safetestset "app" begin
    include("app.jl")
end
@safetestset "reduction visualization" begin
    include("reductionvisualization.jl")
end
@safetestset "aqua" begin
    include("aqua.jl")
end


text:
# =================================================================
# Robertson chemical model
# See
# =================================================================

using ReachabilityAnalysis, Plots

const Tf = 40.0
const normalized_box = symmetric_box(3, Float64)

@taylorize function robertson!(du, u , p, t)
    local α, β, γ = p
    x, y, z = u

    ax = α * x
    yz = y * z
    byz = β * yz
    aux1 = ax - byz
    y2 = y^2
    aux2 = γ * y2

    du[1] = -aux1
    du[2] =  aux1 - aux2
    du[3] =  aux2

    return nothing
end

function robertson()
    # initial state
    X0= (1 .. 1) × (0 .. 0) × (0 .. 0)

    return @ivp(x'= robertson!(x), dim:3, x(0) ∈ X0)
end


text:
getDays(resource,subcal) = (x->x.date[1]).(filter(day ->day.date in subcal,resource.calendar.workdays))

function buildJd(timeslots,subcal)
     timeslots |> @filter(_.dayID in keys(subcal) && !_.booked)|>@groupby(_.resourceID) |>@map((d = key(_),j = unique(map(x->x.dayID,_)))) |> NDSparse
end
function buildI(timeslots)
    timeslots|> @filter(!_.booked) |> @groupby((_.resourceID,_.dayID)) |> @map({d=key(_)[1],j=key(_)[2],i=map(x->x.timeslotID,_)})|> NDSparse
end

function getTlowerbound(timeslots,resourceID,dayID)
    bookedslots = timeslots |> @filter(_.resourceID == resourceID && _.dayID == dayID && _.booked )|> @orderby(_.endTime) |> @take(1)  |>@map(_.endTime) |> collect
    if length(bookedslots) == 1
        value = Dates.value(first(bookedslots))/60000000000
    else
        value = 0
    end
    return value
end

getTSendtime(resource,dayID,slotID) = Dates.value(filter(x->x.date[1]==dayID,resource.calendar.workdays)[1].timeslots[slotID].endTime)/60000000000
function getTSendtime(resource,dayID,slotID,refDate)
    workday = filter(x->x.date[1]==dayID,resource.calendar.workdays)[1]
    datetime = workday.timeslots[slotID].endTime + workday.date[2]
    Dates.value(datetime-refDate)/60000
end
function getTSstarttime(resource,dayID,slotID,refDate)
    workday = filter(x->x.date[1]==dayID,resource.calendar.workdays)[1]
    test = length(workday.timeslots)
    slotID > test ? error("slotID out of range") :
    datetime = workday.timeslots[slotID].startTime + workday.date[2]
    Dates.value(datetime-refDate)/60000
end
getTSstarttime(resource,dayID,slotID) = Dates.value(filter(x->x.date[1]==dayID,resource.calendar.workdays)[1].timeslots[slotID].startTime)/60000000000

function buildVp(visits)
    visits |> @groupby(_.patientID) |> @map((p = key(_),v = map(x->x.intID,_))) |> NDSparse
end
function buildDv(resources,visits)
    visits |> @groupjoin(resources,_.req_type,_.type, (v=Int64(_.intID),d=map(x->x.intID,__))) |> NDSparse
end

function buildDp(resources,visits)
    visits |> @groupjoin(resources,_.req_type,_.type, (v=_.intID,p=_.patientID,d=map(x->x.intID,__))) |> @groupby(_.p) |>@map((p = key(_),d = (collect(Iterators.flatten(_.d)) ))) |>NDSparse
end

"Helper function for grouping patients. Produces a tuple of the patients group"
function patientgroup(visits,Vp,id)
   tempvisits = visits[Vp[id].v]

   (sort(JuliaDB.select(tempvisits,:req_type)),
   Dates.Month(reduce(max,tempvisits,select=:bestord_date)),
   Dates.Month(reduce(min,tempvisits,select=:bestord_date)))
end

"Produces a IndexedTable of patient groups based on the function patientgroup()"
function buildPg(patients,visits,Vp,timeDelta)
    table(patients |> @groupby(patientgroup(visits,Vp,_.intID))|> @map((types=sortvisit(key(_)[1],timeDelta),startmonth=key(_)[2],endmonths=key(_)[3],patients= map(x->x.intID,_))) |> collect)
end

getIndexofPositiveVariables(vars) = filter(k-> value(vars[k]) > 0 ,eachindex(vars))
getvalueofPositiveVariables(vars) = value.((x->vars[x]).(filter(k-> value(vars[k]) > 0 ,eachindex(vars))))
getPositiveVariables(vars) = (x->vars[x]).(filter(k-> value(vars[k]) > 0 ,eachindex(vars)))
getPositiveVariablesAndValues(vars) = (x->(vars[x], value(vars[x]))).(filter(k-> value(vars[k]) > 0 ,eachindex(vars)))

hashSets(subproblem) = hash((subproblem.D_v,subproblem.J_d,subproblem.I,subproblem.Ts,subproblem.Te))


text:
export ScaledEuclidean, ScaledFubiniStudy, FubiniStudyDistances, FubiniStudyDistances!, EllipsoidDistances

"""
EuclideanDistances(data::Array{T,2})

Returns the pairwise euclidean distances of the columns of data, and scales them, such that the maximal pairwise distances is 1.
"""
function EuclideanDistances(data::Array{T,2}) where {T<:Number}
    n = size(data, 2)
    D = zeros(T, n, n)
    D = map(Base.Iterators.product(1:n, 1:n)) do i
        if i[1] < i[2]
            return LinearAlgebra.norm(data[:,i[1]] - data[:,i[2]])
        else
            return 0.0
        end
    end

    return (D+transpose(D))
end
function ScaledEuclidean(data::Array{T,2}) where {T<:Number}
    D = EuclideanDistances(data)
    m = maximum(D)
    return D./m
end

"""
FubiniStudyDistances(data::Array{T,2})

Returns the pairwise Fubini-Study distance of the of the columns of data.
"""
function FubiniStudyDistances(data::Array{T,2}) where {T<:Number}
    m = size(data, 2)
    norms = [LinearAlgebra.norm(data[:,i]) for i in 1:m]

    D = map(Base.Iterators.product(1:m, 1:m)) do i
        if i[1] < i[2]
            p = abs(LinearAlgebra.dot(data[:,i[1]], data[:,i[2]]) / (norms[i[1]] * norms[i[2]]))
            if p > 1
                p = 1.0
            end
            return acos(p)
        else
            return 0.0
        end
    end
    return D+transpose(D)
end
function FubiniStudyDistances!(u::Array{T,2}, data::Array{S,2}) where {T,S<:Number}
    m = size(data, 2)
    norms = [LinearAlgebra.norm(data[:,i]) for i in 1:m]

    for i in Base.Iterators.product(1:m, 1:m)
        if i[1] < i[2]
            p = abs(LinearAlgebra.dot(data[:,i[1]], data[:,i[2]]) / (norms[i[1]] * norms[i[2]]))
            if p > 1
                p = 1.0
            end
            view(u, i[1], i[2]) .= acos(p)
        else
            view(u, i[1], i[2]) .= 0.0
        end
    end
    return u+transpose(u)
end

"""
ScaledFubiniStudy(data::Array{T,2})

Returns the pairwise Fubini-Study distance of the columns of data, and scales them, such that the maximal pairwise distances is 1.
"""
function ScaledFubiniStudy(data::Array{T,2}) where {T<:Number}
    D = FubiniStudyDistances(data)
    m = maximum(D)

    return D ./ m
end

"""
SafeDehomogenization(data::Array{T,2})

Provides affine coordinates for the colums of projective data. The affine patch is safely chosen so that no point lies at infinity.
"""
function SafeDehomogenization(data::Array{T,2}) where {T<:Number}
    n = size(data,1)
    m = size(data,2)
    dots = [data[n,i] for i in 1:m]
    if sum(dots .== 0) == 0
        return hcat([data[1:n-1,i]./dots[i] for i in 1:m]...)
    else
        dots = zeros(T, m)
        u = zeros(T,n)
        check = true
        while check
            u = LinearAlgebra.normalize(rand(n))
            dots = [LinearAlgebra.dot(u, data[:,i]) for i in 1:m]
            if sum(dots .== 0) == 0
                check = false
            end
        end

        S = LinearAlgebra.svd(Diagonal{T}(I, n) - u * u')
        data = hcat([data[:,i]./dots[i] for i in 1:m]...)

        return S.U[:,1:n-1]' * data
    end
end

"""
EllipsoidDistance(data::Array{T,2}, f, λ::Number)

Returns the matrix with distances between points weighted by the radii in direction data[:,i]-data[:,j] of ellipsoids that have radius 1 in tangent direction of {f=0} and radius λ in normal direction of {f=0}.
"""

function EllipsoidDistances(data::Array{T,2}, f::Vector{S}, λ::Number) where {T<:Number, S<:MP.AbstractPolynomialLike}
    n, m = size(data)
    F = FP.System(convert(Vector{FP.Polynomial{Float64}}, f))
    cfg = FP.JacobianConfig(F)
    Js = [FP.jacobian(F, data[:,i], cfg) for i in 1:m]
    Vs = [LinearAlgebra.svd(J, full = true).V for J in Js]
    r = LinearAlgebra.rank(Js[rand(1:m)])
    Σ = LinearAlgebra.diagm(0 => [λ .* ones(r); ones(n-r)])

    Qs = [V * Σ * V' for V in Vs]
    dists = ScaledEuclidean(data)

    D = map(Base.Iterators.product(1:m, 1:m)) do i
     if i[1] < i[2]
          Q1 = Qs[i[1]]
          Q2 = Qs[i[2]]
          u = data[:,i[1]]
          v = data[:,i[2]]

          h = LinearAlgebra.normalize(v-u)

          a1 = LinearAlgebra.transpose(h) * Q1 * h
          a2 = LinearAlgebra.transpose(h) * Q2 * h
          return 2 * dists[i[1],i[2]]/(sqrt(a1[1])+sqrt(a2[1]))
      else
          return 0.0
      end
    end

    for i in findall(isnan.(D))
        D[i[1],i[2]] = 0.0
    end

    return D+transpose(D)
end

function EllipsoidDistances(data::Array{T,2}, f::S, λ::Number) where {T<:Number, S<:MP.AbstractPolynomialLike}
        EllipsoidDistances(data, [f], λ)
end

#
# """
#
# barcode_plot(C::Dict{String,Any},
#     dims::Array{Int64,1},
#     how_many_bars::Array{Int64,1};
#     sorted_by::Symbol=:lower_limit,
#     lw=10,
#     upper_limit = Inf,
#     fontsize = 16)
#
# Plots the barcode associated to the dictionary C that was produced by the eirene() function from the Eirene package. dims is an array determining which dimensions should be plotted. how_many_bars is an array that determines how many bars in each dimension are plotted. sorted_by can be either :lower_limit or :length. It is also possible to pass an array of barcodes B instead of C.
# """
# function barcode_plot(B,
#      dims::Array{Int64,1}, how_many_bars::Array{Int64,1}; sorted_by::Symbol=:lower_limit, lw=10, upper_limit = Inf, fontsize = 16)
#
#         @assert length(dims) == length(how_many_bars) "Number of dimensions (was $(length(dims))) must be the same as the number of values specifying the number of bars that should be displayed for each dimension  (was $(length(how_many_bars)))."
#
#         # colors = Colors.distinguishable_colors(length(dims)+1,[RGB(1,1,1)])[2:end]
#         cols = Colors.colormap("Blues", mid = 0.5)
#         range = Int.(round.(collect(linspace(50,100,length(dims)+1))))
#         colors = map(r -> cols[r], range)
#
#         if upper_limit == Inf
#             upper_limit = 2 * maximum([maximum(b[b.< Inf]) for b in B])
#         end
#
#         traces = map(1:length(B)) do j
#              b = B[j]
#              # lengths = b[:,2]-b[:,1]
#              # b = b[lengths .> tol[j],:]
#              if size(b,1) == 0
#                  return [PlotlyJS.scatter(;x=[0,0], y=[0,0], mode="lines",  line_width = lw, line_color = colors[j], name = "dimension $(dims[j])")]
#              end
#
#              l = minimum([size(b,1), how_many_bars[j]])
#              s = sortperm(b[:,2]-b[:,1])
#              b = b[s[end-l+1:end],:]
#
#              i = find(x->x==Inf, b[:,2])
#              b[i,2] .= upper_limit
#
#              if sorted_by == :length
#              elseif sorted_by == :lower_limit
#                 s = sortperm(b[:,1])
#                 b = b[s,:]
#              else
#                 println("The second argument must be either :length or :lower_limit.")
#                 return 0
#              end
#              return [PlotlyJS.scatter(;x=b[1,:], y=[1,1], mode="lines",  line_width = lw, line_color = colors[j], name = "Dimension $(dims[j])");
#              [PlotlyJS.scatter(;x=b[i,:], y=[i,i], mode="lines",  line_width = lw, line_color = colors[j], showlegend=false) for i in 2:l]]
#
#          end
#
#          for i = 2:length(traces)
#              for j = 1:length(traces[i])
#                  traces[i][j][:y] = traces[i][j][:y] + traces[i-1][end][:y] + 1
#              end
#          end
#          traces = vcat(traces...)
#
#
#          x = maximum(vcat([t[:x] for t in traces]...))
#          y = maximum(vcat([t[:y] for t in traces]...))
#
#          layout = PlotlyJS.Layout(;
#          xaxis = attr(range = [-.001, upper_limit+0.001], showgrid=false, zeroline =false, title = "ϵ", titlefont_size=fontsize, tickfont_size=fontsize),
#          yaxis = attr(range = [0,y+1.1], showgrid=false, ticks = false, titlefont_size=fontsize, tickfont_size=fontsize),
#          legend=attr(font_size=fontsize)
#          )
#          PlotlyJS.Plot(traces[end:-1:1], layout)
# end
# function barcode_plot(C::Dict{String,Any}, dims::Array{Int64,1}, how_many_bars::Array{Int64,1}; sorted_by::Symbol=:lower_limit, lw=10, upper_limit = Inf, fontsize = 16)
#     B = [barcode(C, dim = d) for d in dims]
#     barcode_plot(B, dims, how_many_bars, sorted_by=sorted_by, lw=lw, upper_limit = upper_limit, fontsize = fontsize)
# end
# function barcode_plot(C::Dict{String,Any}, dims::Array{Int64,1}; sorted_by::Symbol=:lower_limit, lw=10, upper_limit = Inf, fontsize = 16)
#     l = size(C["symmat"],1)
#     B = [barcode(C, dim = d) for d in dims]
#     barcode_plot(B, dims, l .* ones(Int64,length(dims)), sorted_by=sorted_by, lw=lw, upper_limit = upper_limit, fontsize = fontsize)
# end