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

import DelimitedFiles

# ==================== GATHERING METHODS ====================

function gather_data_from_files(
	@nospecialize(filenames), @nospecialize(extractors)
)
	file_contents = String[]
	sizehint!(file_contents, length(filenames))
	for filename in filenames
		push!(file_contents, read(filename, String))
	end
	processed_data = Any[]
	sizehint!(processed_data, length(extractors))
	for e in extractors
		push!(processed_data, map(e, file_contents))
	end
	return processed_data
end

function gather_data_from_folder(
	@nospecialize(folder_name), @nospecialize(extractors)
)
	filenames = readdir(folder_name; join = true, sort = false)
	return gather_data_from_files(filenames, extractors)
end

function _data2csv(
	@nospecialize(data), @nospecialize(delim), @nospecialize(column_names)
)
	iob = IOBuffer()
	if column_names !== missing
		println(iob, join(column_names, delim))
	end
	DelimitedFiles.writedlm(iob, zip(data...), delim)
	csv_str = read(seekstart(iob), String)
	return csv_str
end

function gather_csv_from_folders(
	@nospecialize(folders_name), @nospecialize(extractors);
	delim = ';', column_names = missing
) :: String
	if !ismissing(column_names)
		@assert length(column_names) == length(extractors)
	end
	data = Any[]
	for folder_name in folders_name
		push!(data, gather_data_from_folder(folder_name, extractors))
	end
	#if isone(length(data))
	#	data = only(data)
	#else
		data = map(vcat, data...)
	#end
	return _data2csv(data, delim, column_names)
end

function gather_csv_from_folder(
	@nospecialize(folder_name), @nospecialize(extractors);
	delim = ';', column_names = missing
) :: String
	return gather_csv_from_folders(
		[folder_name], extractors; delim = delim, column_names = column_names
	)
end

# ==================== MY_PARSE METHODS ====================

# An extensible parsing method.
function my_parse(::Type{String}, str :: AbstractString) :: String
	m = match(r"^\"(.*)\"$", str)
	if isnothing(m)
		return String(str)
	else
		return String(m.captures[1])
	end
end

function my_parse(::Type{N}, str :: AbstractString) :: N where {N <: Number}
	return parse(N, str)
end

# ==================== EXTRACTOR METHODS ====================

struct NoDefault{T} end

struct RegexExtractor{T}
	regex :: Regex
	default :: Union{T, NoDefault{T}}
end

function (re::RegexExtractor{T})(data :: AbstractString) :: T where {T}
	m = match(re.regex, data)
	isnothing(m) || return my_parse(T, only(m.captures))
	if re.default === NoDefault{T}()
		error("The regex $(re.regex) has found no matches. No default available.")
	end
	return re.default
end

function key_equals_extractor(
	key :: AbstractString, default :: Union{T, NoDefault{T}}
) where {T}
	return RegexExtractor{T}(
		Regex("^$key = (.*)\$", "m"), default
	)
end

function p_args_key_extractor(
	key :: AbstractString, default :: Union{T, NoDefault{T}}
) where {T}
	return RegexExtractor{T}(
		Regex("^p_args = .*\"$key\" => ([^,)]*)", "m"), default
	)
end

struct BooleanExtractor
	regex :: Regex
	match_answer :: Bool
end

function (be::BooleanExtractor)(data :: AbstractString)
	m = match(be.regex, data)
	isnothing(m) && return !be.match_answer
	return be.match_answer
end

function matches(r :: Regex)
	return BooleanExtractor(r, true)
end

function does_not_match(r :: Regex)
	return BooleanExtractor(r, false)
end

@enum RNStat Objective=1 Iterations=2 Time=3
Base.getindex(indexable :: AbstractArray, s :: RNStat) = indexable[Int(s)]
Base.getindex(indexable :: Tuple, s :: RNStat) = indexable[Int(s)]

struct RootNodeStatsExtractor
	marker :: String
	stat :: RNStat
end

# Return the objective value, number of iterations, and the time spent
# (in seconds) in the first Gurobi log node relaxation line it finds after
# the `marker` string (the `maker` should be a full line).
function (rnse::RootNodeStatsExtractor)(
	log :: AbstractString
) :: Union{Int, Float64}
	marker = "MARK_" * rnse.marker
	# First, let us restrict the log to the part after our marker and before the
	# next marker (if there is one, otherwise until the end of the input).
	all_markers_positions = findall(r"^MARK_.+$"m, log)
	all_markers = getindex.(log, all_markers_positions)
	our_marker_idx = findfirst(isequal(marker), all_markers)
	isnothing(our_marker_idx) && return (NaN, -1, NaN)[rnse.stat]
	first_byte = first(all_markers_positions[our_marker_idx])
	if our_marker_idx == length(all_markers)
		# If it is the last marker, just go from there to the end of the string.
		final_byte = length(log)
	else
		# If it is not the last marker, get from there to the start of the next.
		# The final byte is always an 'M' character in this case.
		final_byte = first(all_markers_positions[our_marker_idx + 1])
	end
	log = log[first_byte:final_byte]
	# The presolve can remove all rows and columns in some cases, and then
	# the extraction below will fail. So we check this first, and return zero
	# values in this case (that is the most sensible value).
	presolve_regex = Regex("^Presolve: All rows and columns removed\$", "m")
	if !isnothing(match(presolve_regex, log))
		return (0.0, 0, 0.0)[rnse.stat]
	end
	# If barrier is used, the barrier log line appears first and has the most
	# accurate time. If it isn't the root relaxation line is the only one
	# present and has the correct time. Unfortunately the order of the captures
	# in this segment does not match the RNStat enum.
	barrier_log = "^Barrier solved model in (\\d+) iterations and (\\S+) " *
		"seconds\nOptimal objective (\\S+)"
	# The order of captures in the line below match the enum RNStat.
	root_rel = "^Root relaxation: (?:objective (\\S+)|cutoff), (\\d+) " *
		"iterations, (\\S+) seconds"
	either_line = "(?:$barrier_log|$root_rel)"
	r = Regex(either_line, "m")
	m = match(r, log)
	if nothing === m
		return (NaN, -1, NaN)[rnse.stat]
	else
		@assert 6 == length(m.captures)
		# If there is no barrie line.
		if isnothing(match(r"^Barrier solved model"m, m.match))
			# A capture has index i if it is the i-esim capture to appear in the
			# string, even if previous captures cannot coexist with the current
			# one (i.e., one is before a '|' inside a group and the other after).
			# This way, all non-barrier captures appear after the barrier ones,
			# even if there were no barrier section (the first three captures
			# are nothing in this case).
			@assert m.captures[[1, 2, 3]] == [nothing, nothing, nothing]
			idx = 3 + Int(rnse.stat) # The enumeration has the right order.
			# The root node relaxation can present the string 'cutoff' instead
			# of the objective value, not sure, but it seems to happen when the
			# warm-start has the same value as the relaxation rounded down. In
			# this case, the script processing the data needs to search the
			# objective value in the restricted_LP key. We will not do this
			# here.
			if rnse.stat == Objective && isnothing(m.captures[idx])
				@assert !isnothing(match(r"^Root relaxation: cutoff,"m, m.match))
				return NaN
			end
			return parse((Float64, Int, Float64)[rnse.stat], m.captures[idx])
		else # If there is a barrier line.
			@assert m.captures[[4, 5, 6]] == [nothing, nothing, nothing]
			# The enumeration has not the right order in this line.
			idx :: Int = rnse.stat == Iterations ? 1 : rnse.stat == Time ? 2 : 3
			return parse((Int, Float64, Float64)[idx], m.captures[idx])
		end
	end
end

# ==================== FUNCTION CALLS ====================

# Data extraction for the experiment over the Velasco and Uchoa instances.
#=
csv = gather_csv_from_folders(
	"./finished_experiments/" .* [
		"vu_ixion_2020-07-28T05:01:05/",
		"vu_leviathan_2020-07-28T04:59:52/",
		"vu_ramuh_2020-07-28T05:02:27/",
		"vu_yojimbo_2020-07-28T05:00:34/",
		"vu_odin_2020-07-28T05:58:36/"
	],
	[
		key_equals_extractor("instfname", NoDefault{String}()),
		p_args_key_extractor("PPG2KP-pricing", NoDefault{String}()),
		key_equals_extractor("qt_pevars_priced_restricted", -1),
		key_equals_extractor("qt_cmvars_priced_restricted", -1),
		key_equals_extractor("qt_plates_priced_restricted", -1),
		key_equals_extractor("qt_pevars_after_purge", -1),
		key_equals_extractor("qt_cmvars_after_purge", -1),
		key_equals_extractor("qt_plates_after_purge", -1),
		key_equals_extractor("heuristic_lb", NaN),
		key_equals_extractor("restricted_obj_value", NaN),
		key_equals_extractor("restricted_obj_bound", NaN),
		key_equals_extractor("restricted_lp_stop_reason", "NOT_REACHED"),
		key_equals_extractor("restricted_stop_reason", "NOT_REACHED"),
		key_equals_extractor("solved_priced_restricted_model", false),
		key_equals_extractor("enumeration_time", NaN),
		key_equals_extractor("heuristic_lb_time", NaN),
		key_equals_extractor("restricted_pricing_time", NaN),
		key_equals_extractor("iterative_pricing_time", NaN),
		key_equals_extractor("final_pricing_time", NaN),
		key_equals_extractor("pricing_time", NaN),
		RootNodeStatsExtractor("FINAL_GENERIC_SOLVE", Time),
		key_equals_extractor("finished_model_solve", NaN),
		key_equals_extractor("obj_value", NaN),
		key_equals_extractor("solution_profit", NaN),
		key_equals_extractor("obj_bound", NaN),
		key_equals_extractor("total_instance_time", NaN),
		key_equals_extractor("build_and_solve_time", NaN),
		key_equals_extractor("run_total_time", NoDefault{Float64}()),
		key_equals_extractor("run_ended_by_exception", NoDefault{Bool}()),
		matches(r"TimeoutError"),
		key_equals_extractor("build_stop_reason", "NOT_REACHED"),
		key_equals_extractor("this_data_file", NoDefault{String}())
	];
	column_names = [
		"instance_name",
		"pricing_method",
		"qt_pevars_priced_restricted",
		"qt_cmvars_priced_restricted",
		"qt_plates_priced_restricted",
		"qt_pevars_after_purge",
		"qt_cmvars_after_purge",
		"qt_plates_after_purge",
		"heuristic_lb",
		"restricted_obj_value",
		"restricted_obj_bound",
		"restricted_lp_stop_reason",
		"restricted_stop_reason",
		"solved_priced_restricted_model",
		"enumeration_time",
		"heuristic_lb_time",
		"restricted_pricing_time",
		"iterative_pricing_time",
		"final_pricing_time",
		"total_pricing_time",
		"final_root_relaxation_time",
		"final_solving_time",
		"obj_value",
		"solution_profit",
		"obj_bound",
		"total_instance_time",
		"build_and_solve_time",
		"run_total_time",
		"run_ended_by_exception",
		"had_timeout",
		"build_stop_reason",
		"datafile"
	]
)
print(csv)
=#

# Data extraction for the experiment comparing the revised model with
# the our reimplementation of the original model.
#=
csv = gather_csv_from_folders(
	"./finished_experiments/comparison_ramuh/" .* [
		"2020-08-03T14:56:36/",
		"2020-08-04T13:05:52/",
		"2020-08-07T14:27:24/",
		"2020-08-10T13:19:35/",
		"2020-08-11T13:25:03/",
	],
	[
		key_equals_extractor("instfname", NoDefault{String}()),
		p_args_key_extractor("PPG2KP-pricing", NoDefault{String}()),
		p_args_key_extractor("PPG2KP-round2disc", NoDefault{String}()),
		p_args_key_extractor("PPG2KP-MIP-start", NoDefault{String}()),
		p_args_key_extractor(
			"PPG2KP-do-not-purge-unreachable", NoDefault{String}()
		),
		p_args_key_extractor(
			"PPG2KP-faithful2furini2016", NoDefault{String}()
		),
		key_equals_extractor("build_and_solve_time", NaN),
		key_equals_extractor("total_instance_time", NaN),
		key_equals_extractor("final_pricing_time", NaN),
		key_equals_extractor("iterative_pricing_time", NaN),
		key_equals_extractor("restricted_pricing_time", NaN),
		key_equals_extractor("heuristic_lb_time", NaN),
		key_equals_extractor("enumeration_time", NaN),
		key_equals_extractor("pricing_time", NaN),
		RootNodeStatsExtractor("FINAL_GENERIC_SOLVE", Time),
		key_equals_extractor("finished_model_solve", NaN),
		key_equals_extractor("qt_pevars_after_preprocess", -1),
		key_equals_extractor("qt_cmvars_after_preprocess", -1),
		key_equals_extractor("qt_plates_after_preprocess", -1),
		key_equals_extractor("length_pe_after_pricing", -1),
		key_equals_extractor("length_cm_after_pricing", -1),
		key_equals_extractor("length_pc_after_pricing", -1),
		key_equals_extractor("qt_pevars_after_purge", -1),
		key_equals_extractor("qt_cmvars_after_purge", -1),
		key_equals_extractor("qt_plates_after_purge", -1),
		key_equals_extractor("run_total_time", NoDefault{Float64}()),
		key_equals_extractor("solution_profit", NaN),
		# The default of 'run_ended_by_exception' is that the code has
		# ended by exception. The reasoning is that some exceptions are not
		# capturable, and so the absence of this field indicates an exception
		# happened.
		key_equals_extractor("run_ended_by_exception", true),
		matches(r"TimeoutError"),
		key_equals_extractor("build_stop_reason", "NOT_REACHED"),
		key_equals_extractor("this_data_file", NoDefault{String}())
	];
	column_names = [
		"instance_name",
		"pricing_method",
		"round2disc",
		"mip_start",
		"purge_disabled",
		"faithful",
		"build_and_solve_time",
		"total_instance_time",
		"final_pricing_time",
		"iterated_pricing_time",
		"restricted_pricing_time",
		"heuristic_lb_time",
		"enumeration_time",
		"total_pricing_time",
		"final_root_relaxation_time",
		"final_solving_time",
		"qt_pevars_after_preprocess",
		"qt_cmvars_after_preprocess",
		"qt_plates_after_preprocess",
		"qt_pevars_after_pricing",
		"qt_cmvars_after_pricing",
		"qt_plates_after_pricing",
		"qt_pevars_after_purge",
		"qt_cmvars_after_purge",
		"qt_plates_after_purge",
		"run_total_time",
		"solution_profit",
		"run_ended_by_exception",
		"had_timeout",
		"build_stop_reason",
		"datafile"
	]
)
print(csv)
=#

# Data extraction for the experiment related to barrier vs dual simplex
# and their effects in the furini pricing.
#=
csv = gather_csv_from_folder(
	"./finished_experiments/LP_method_2020-07-31T19:08:28/",
	[
		key_equals_extractor("instfname", NoDefault{String}()),
		p_args_key_extractor("Gurobi-LP-method", NoDefault{Int}()),
		p_args_key_extractor(
			"PPG2KP-Gurobi-LP-method-inside-furini-pricing", NoDefault{Int}()
		),
		p_args_key_extractor("PPG2KP-pricing", NoDefault{String}()),
		RootNodeStatsExtractor("FURINI_PRICING_RESTRICTED_MIP_SOLVE", Time),
		RootNodeStatsExtractor("FINAL_GENERIC_SOLVE", Time),
		key_equals_extractor("build_and_solve_time", NaN),
		key_equals_extractor("total_instance_time", NaN),
		key_equals_extractor("final_pricing_time", NaN),
		key_equals_extractor("iterative_pricing_time", NaN),
		key_equals_extractor("restricted_pricing_time", NaN),
		key_equals_extractor("pricing_time", NaN),
		key_equals_extractor("solution_profit", NaN),
		matches(r"TimeoutError"),
		# It is kinda terrible to depend on this, but the type information of the
		# printed backtrace is the most reliable indicator that the run ended by
		# exception. New runs will have a better unequivocal indicator.
		does_not_match(r"StackTraces"),
		key_equals_extractor("this_data_file", NoDefault{String}())
	];
	column_names = [
		"instance_name",
		"lp_method",
		"lp_method_switch",
		"pricing_method",
		"restricted_root_time",
		"final_root_time",
		"build_and_solve_time",
		"total_instance_time",
		"final_pricing_time",
		"iterated_pricing_time",
		"restricted_pricing_time",
		"total_pricing_time",
		"solution_profit",
		"had_timeout",
		"finished",
		"datafile"
	]
)
print(csv)
=#

# TODO: to run this we need to solve the problem with the two gcut12
# runs that had memory problems.
# Faithful reimplementation data.
#=
csv = gather_csv_from_folders(
	[
		"./finished_experiments/faithful_2020-07-31T21:24:07/",
		"./finished_experiments/faithful_2020-08-03T18:32:40/",
		"./finished_experiments/faithful_2020-08-04T12:48:50/",
		"./finished_experiments/faithful_2020-08-04T18:41:05/",
	],
	[
		# Primary key: these four identify a set of parameters and, as there
		# was no repetition, a row of the CSV.
		key_equals_extractor("instfname", NoDefault{String}()),
		p_args_key_extractor("PPG2KP-pricing", NoDefault{String}()),
		p_args_key_extractor("PPG2KP-no-redundant-cut", NoDefault{Bool}()),
		p_args_key_extractor("PPG2KP-no-cut-position", NoDefault{Bool}()),

		# The model stats at multiple stats.
		key_equals_extractor("qt_pevars_after_preprocess", -1),
		key_equals_extractor("qt_cmvars_after_preprocess", -1),
		key_equals_extractor("qt_plates_after_preprocess", -1),

		key_equals_extractor("qt_pevars_restricted", -1),
		key_equals_extractor("qt_cmvars_restricted", -1),
		key_equals_extractor("qt_plates_restricted", -1),

		key_equals_extractor("qt_pevars_priced_restricted", -1),
		key_equals_extractor("qt_cmvars_priced_restricted", -1),
		key_equals_extractor("qt_plates_priced_restricted", -1),

		key_equals_extractor("qt_pevars_before_iterated", -1),
		key_equals_extractor("qt_cmvars_before_iterated", -1),
		key_equals_extractor("qt_plates_before_iterated", -1),

		key_equals_extractor("qt_pevars_after_iterated", -1),
		key_equals_extractor("qt_cmvars_after_iterated", -1),
		key_equals_extractor("qt_plates_after_iterated", -1),

		key_equals_extractor("qt_pevars_after_final", -1),
		key_equals_extractor("qt_cmvars_after_final", -1),
		key_equals_extractor("qt_plates_after_final", -1),

		key_equals_extractor("length_pe_after_pricing", -1),
		key_equals_extractor("length_cm_after_pricing", -1),
		key_equals_extractor("length_pc_after_pricing", -1),

		key_equals_extractor("qt_pevars_after_purge", -1),
		key_equals_extractor("qt_cmvars_after_purge", -1),
		key_equals_extractor("qt_plates_after_purge", -1),

		key_equals_extractor("build_stop_reason", "NOT_REACHED"),
		matches(r"TimeoutError"),
		# The default of 'run_ended_by_exception' is that the code has
		# ended by exception. The reasoning is that some exceptions are not
		# capturable, and so the absence of this field indicates an exception
		# happened.
		key_equals_extractor("run_ended_by_exception", true),
		key_equals_extractor("this_data_file", NoDefault{String}()),
	];
	column_names = [
		# primary key
		"instance_name",
		"pricing_method",
		"disabled_redundant_cut",
		"disabled_cut_position",
		# The model stats at multiple stats.
		"qt_pevars_after_preprocess",
		"qt_cmvars_after_preprocess",
		"qt_plates_after_preprocess",

		"qt_pevars_restricted",
		"qt_cmvars_restricted",
		"qt_plates_restricted",

		"qt_pevars_priced_restricted",
		"qt_cmvars_priced_restricted",
		"qt_plates_priced_restricted",

		"qt_pevars_before_iterated",
		"qt_cmvars_before_iterated",
		"qt_plates_before_iterated",

		"qt_pevars_after_iterated",
		"qt_cmvars_after_iterated",
		"qt_plates_after_iterated",

		"qt_pevars_after_final",
		"qt_cmvars_after_final",
		"qt_plates_after_final",

		# The names of the three stats below are normalized.
		"qt_pevars_after_pricing", # "length_pe_after_pricing",
		"qt_cmvars_after_pricing", # "length_cm_after_pricing",
		"qt_plates_after_pricing", # "length_pc_after_pricing",

		"qt_pevars_after_purge",
		"qt_cmvars_after_purge",
		"qt_plates_after_purge",
		# extra info
		"build_stop_reason",
		"had_timeout",
		"run_ended_by_exception",
		"datafile",
	]
)
print(csv)
=#

# Data extraction for the LAGOS 2021.
path = "raw_data/" # inside latex/lagos2021/notebooks
csv = gather_csv_from_folders(
	path .* vcat(
		"g2csp/" .* [
			"dataset_c_leviathan_2020-11-16T19:13:37",
			"dataset_c_leviathan_2020-11-18T20:28:18",
			"dataset_c_odin_2020-11-16T19:01:25",
			"dataset_c_yojimbo_2020-11-17T03:16:22",
			"gcuts_yojimbo_2020-11-17T16:09:01",
			"odin_extra_seed_2020-11-23T21:32:16",
			"bahamut_extra_seeds_2020-11-23T21:29:54",
		],
		"g2kp/" .* [
			"dataset_c_bahamut_2020-11-17T03:18:25",
			"dataset_c_bahamut_2020-11-18T18:49:04",
			"gcuts_bahamut_2020-11-19T13:59:09",
		]
	),
	[
		key_equals_extractor("instfname", NoDefault{String}()),
		p_args_key_extractor("PPG2KP-round2disc", NoDefault{String}()),
		p_args_key_extractor(
			"PPG2KP-faithful2furini2016", NoDefault{String}()
		),
		p_args_key_extractor(
			"Gurobi-threads", NoDefault{Int}()
		),
		p_args_key_extractor(
			"Gurobi-seed", NoDefault{Int}()
		),
		key_equals_extractor("build_and_solve_time", NaN),
		key_equals_extractor("total_instance_time", NaN),
		key_equals_extractor("enumeration_time", NaN),
		matches(r"Barrier solve interrupted - model solved by another algorithm"),
		matches(r"Solved with primal simplex"),
		matches(r"Solved with dual simplex"),
		RootNodeStatsExtractor("FINAL_GENERIC_SOLVE", Time),
		RootNodeStatsExtractor("FINAL_GENERIC_SOLVE", Objective),
		key_equals_extractor("finished_model_solve", NaN),
		key_equals_extractor("qt_pevars_after_preprocess", -1),
		key_equals_extractor("qt_cmvars_after_preprocess", -1),
		key_equals_extractor("qt_plates_after_preprocess", -1),
		key_equals_extractor("qt_pevars_after_purge", -1),
		key_equals_extractor("qt_cmvars_after_purge", -1),
		key_equals_extractor("qt_plates_after_purge", -1),
		key_equals_extractor("run_total_time", NoDefault{Float64}()),
		key_equals_extractor("solution_profit", -1), # if G2KP
		key_equals_extractor("solution_value", -1), # if G2CSP
		# The default of 'run_ended_by_exception' is that the code has
		# ended by exception. The reasoning is that some exceptions are not
		# capturable and, when this happens, the code that writes this field
		# does not run, so the absence of this field indicates an exception
		# happened.
		key_equals_extractor("run_ended_by_exception", true),
		matches(r"TimeoutError"),
		key_equals_extractor("build_stop_reason", "NOT_REACHED"),
		key_equals_extractor("this_data_file", NoDefault{String}())
	];
	column_names = [
		"instance_name",
		"round2disc",
		"faithful",
		"threads",
		"seed",
		"build_and_solve_time",
		"total_instance_time",
		"enumeration_time",
		"rnr_solved_by_simplex",
		"rnr_solved_by_primal",
		"rnr_solved_by_dual",
		"final_root_relaxation_time",
		"final_root_relaxation_value",
		"final_solving_time",
		"qt_pevars_after_preprocess",
		"qt_cmvars_after_preprocess",
		"qt_plates_after_preprocess",
		"qt_pevars_after_purge",
		"qt_cmvars_after_purge",
		"qt_plates_after_purge",
		"run_total_time",
		"solution_profit",
		"solution_value",
		"run_ended_by_exception",
		"had_timeout",
		"build_stop_reason",
		"datafile"
	]
)
print(csv)


text:
# ## CDPlayer

# ## Model

module cdplayer #jl

using ReachabilityAnalysis, MAT
using ReachabilityModels: @relpath

file = matopen(@relpath "cdplayer.mat")

## system matrix
A = read(file, "A")

## input matrix
B = read(file, "B")

## state domain
X = Universe(120)

## input domain
U = BallInf([0.0, 0.0], 1.0)

function model(X0)
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)
    return IVP(S, X0)
end

end  #jl

# ## References

ext: jl

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# Copyright (c) 2019 Mathieu Besançon, Oscar Dowson, and contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module HiGHS

import HiGHS_jll
import MathOptInterface
import SparseArrays

function __init__()
    global libhighs = HiGHS_jll.libhighs
    return
end

include("gen/libhighs.jl")
include("MOI_wrapper.jl")

# HiGHS exports all `Highs_xxx` symbols. If you don't want all of these symbols
# in your environment, then use `import HiGHS` instead of `using HiGHS`.

for sym in names(@__MODULE__, all = true)
    if startswith("$sym", "Highs_") || startswith("$sym", "kHighs")
        @eval export $sym
    end
end

end # module HiGHS


text:
function test_scope()
    local eigen_v;
    for i = 1:5
        eigen_v = 1;
    end
    return eigen_v
end

function test_hello()
    print("hello")
end

function test_emacs()
    print("hello")
end


text:
function extract_nodes(file, ::Type{TF}=Float64, ::Type{TI}=Int) where {TF,TI}
    line = readline(file)

    pattern = r"(\d+)\s*,\s*(-?\d+\.?\d*(e[-\+]?\d*)?)\s*,\s*(-?\d+\.?\d*(e[-\+]?\d*)?)\s*(,\s*(-?\d+\.?\d*(e[-\+]?\d*)?))?"
    m = match(pattern, line)

    first_node_idx = parse(TI, m[1])
    first_node_idx == TI(1) || throw("First node index is not 1.")

    if m[6] isa Nothing
        node_coords = [(parse(TF, m[2]), parse(TF, m[4]))]
        nextline = _extract_nodes!(node_coords, file, TI(1))
    else
        node_coords = [(parse(TF, m[2]), parse(TF, m[4]), parse(TF, m[7]))]
        nextline = _extract_nodes!(node_coords, file, TI(1))
    end
    return node_coords, nextline
end

function _extract_nodes!(
    node_coords::AbstractVector{NTuple{dim,TF}}, file, prev_node_idx::TI
) where {dim,TF,TI}
    if dim === 2
        pattern = r"(\d+)\s*,\s*(-?\d+\.?\d*(e[-\+]?\d*)?)\s*,\s*(-?\d+\.?\d*(e[-\+]?\d*)?)"
    elseif dim === 3
        pattern = r"(\d+)\s*,\s*(-?\d+\.?\d*(e[-\+]?\d*)?)\s*,\s*(-?\d+\.?\d*(e[-\+]?\d*)?)\s*,\s*(-?\d+\.?\d*(e[-\+]?\d*)?)"
    else
        error("Dimension is not supported.")
    end

    line = readline(file)
    m = match(stopping_pattern, line)
    while m isa Nothing
        m = match(pattern, line)
        if m != nothing
            node_idx = parse(Int, m[1])
            node_idx == prev_node_idx + TI(1) || throw("Node indices are not consecutive.")
            if dim === 2
                push!(node_coords, (parse(TF, m[2]), parse(TF, m[4])))
            else
                push!(node_coords, (parse(TF, m[2]), parse(TF, m[4]), parse(TF, m[6])))
            end
            prev_node_idx = node_idx
        end
        line = readline(file)
        m = match(stopping_pattern, line)
    end
    return line
end


text:
# goes through an expression and replace all QuoteNode's that are not in an expression with
# head = ref by an expression that ref's the i'th element
function ref_quotenodes!(e::Expr, headtype::Symbol = :top)
    if (e.head != :ref)
        for aind = 1:length(e.args)
            a = e.args[aind]
            if MacroTools.isexpr(a)
                ref_quotenodes!(a, e.head)
            elseif (typeof(a) == QuoteNode) && (headtype != :ref)
                # replace this
                e.args[aind] = Expr(:ref, a, :_n)
            end
        end
    end
    e
end

# recursively goes through an expression and replaces all QuoteNode's by s[#], while populating
# qn_vec with exactly these QuoteNodes
function replace_QuoteNodes!(e::Expr, s::Symbol, qn_vec::Vector{QuoteNode})
    for aind = 1:length(e.args)
        a = e.args[aind]
        if MacroTools.isexpr(a)
            replace_QuoteNodes!(a, s, qn_vec)
        elseif (typeof(a) == QuoteNode)
            # replace this
            idx = findfirst(isequal(a), qn_vec)
            if !isnothing(idx)
                # already exists in qn_vec
                e.args[aind] = Expr(:ref, s, idx)
            else
                push!(qn_vec, e.args[aind])
                e.args[aind] = Expr(:ref, s, length(qn_vec) )
            end
        end
    end
end

# recursively go through the expression and replace arg[#] by arg[#][index]
function push_index!(e::Expr, index::Symbol, s::Symbol)
    if (e.head == :ref) && (length(e.args)>1) && (e.args[1] == s)
        e.args[1] = copy(e)
        e.args[2] = index
    else
        # go through arguments and if they are expressions, go deeper
        for aind = 1:length(e.args)
            a = e.args[aind]
            if MacroTools.isexpr(a)
                push_index!(a, index, s)
            end
        end
    end
end

# removes the colon at the start of a string, if present
function stripcolon(s::AbstractString)
    return (length(s)>1 && s[1] == ':') ? s[2:end] : s
end

# replace all Expr that take the following form by an expression that evaluates to `missing` if the index is negative
# head: Symbol ref
# args: args: Array{Any}((2,))
#   1: Symbol y
#   2: Expr
function replace_invalid_indices!(ex::Expr)
    if (ex.head == :ref) && isa(ex.args[1], QuoteNode) && MacroTools.isexpr(ex.args[2])
        r = deepcopy(ex)
        myex = :( ($(ex.args[2]) < 1) ? missing : $(r) )
        ex.head = deepcopy(myex.head)
        ex.args = deepcopy(myex.args)
    else
        # keep going through the tree
        for a in ex.args
            MacroTools.isexpr(a) && replace_invalid_indices!(a)
        end
    end
    return ex
end

# returns `true` if `keys` uniquely identify rows in `df`, otherwise `false`
function unique_obs(df::DataFrame, keys::Vector{Symbol})
    return allunique(Tables.namedtupleiterator(df[!,keys]))
end

# expand the argument x to the length of the df if it's not already a vector
# first generic version that supports size(_,1)
function helper_expand(df, x)
    (ismissing(x) || size(x,1) == 1) ? repeat([x],size(df,1)) : x
end
# ... or to a length of l::Int64
function helper_expand(l::Int64, x)
    (ismissing(x) || size(x,1) == 1) ? repeat([x],l) : x
end

function assign_helper!(var, indices, assignee::Vector{T}) where T<:Any
    var[indices] = assignee
end
function assign_helper!(var, indices, assignee::T) where T<:Any
    var[indices] .= assignee
end
# argument with vector
function assign_helper_gen(indices, assignee::Vector{T}) where T<:Any
    var = missings(T,length(indices))
    var[indices] = assignee
    return var
end
# argument with scalar
function assign_helper_gen(indices, assignee::T) where T<:Any
    var = missings(T,length(indices))
    var[indices] .= assignee
    return var
end
# argument with vector
function assign_helper_rep!(var, indices, assignee::Vector{T}) where T<:Any
    var[indices] = assignee
    return var
end
# argument with scalar
function assign_helper_rep!(var,indices, assignee::T) where T<:Any
    var[indices] .= assignee
    return var
end

# Some helper macros

# helper macro to make sure filter is valid
macro assert_filter(t::Symbol, filter::Expr)
    esc(
        quote
            local x = Douglass.@with($t, $filter)
            (typeof(x) <: BitArray{1}) || (typeof(x) <: Vector{Union{Missing,Bool}}) || error("filter is not a valid boolean vector.")
            true
        end
    )
end

# helper macro to make sure that expression evaluates to a Vector{Symbol}
macro assert_varlist(t::Symbol, varlist::Expr)
    esc(
        quote
            typeof($varlist) == Vector{Symbol} || error("Argument must be evaluating to a Vector{Symbol}. Type is $(typeof($varlist))")
            true
        end
    )
end

# asserts that varlist::Expr evaluates to a Vector{Symbol} and checks that all Symbols are column names in t.
macro assert_vars_present(t::Symbol, varlist::Expr)
    esc(
        quote
            # check that it's a varlist
            Douglass.@assert_varlist($t, $varlist)
            # check that they're present
            for v in $varlist
                (v ∈ propertynames($t)) || error("$(v) not a column name in the active DataFrame")
            end
            true
        end
    )
end

# Checks that all Symbols are column names in t.
macro assert_vars_present(t::Symbol, varlist::Vector{Symbol})
    esc(
        quote
            # check that they're present
            ($varlist ⊆ propertynames($t)) || error("$($varlist) is not a subset of the columns in the active DataFrame")
            true
        end
    )
end


text:
using TerminalUI
using Reactive
using VT100
using Compose

import TerminalUI: Border, TightCentering, SimpleSlider,
    DoubleBufferedTerminalScreen, draw, CheckBox, WidgetStack,
    Gauge, RowLayout, MultiLineInput, IOBufferView, ScrollableChain,
    ImageWidget
import Base.Terminals: CSI

checkboxes = [
    CheckBox("Pineapple"),
    CheckBox("Onion"),
    CheckBox("Cheese"),
    CheckBox("Tomato Saunce")
]

tty = Base.Terminals.TTYTerminal("xterm",STDIN,STDOUT,STDERR)

function run_repl(repl)
    repl_channel = Channel(1)
    response_channel = Channel(1)
    backend = Base.REPL.start_repl_backend(repl_channel, response_channel)
    repl.mi.backendref = Base.REPL.REPLBackendRef(repl_channel,response_channel)
    backend
end

sl = SimpleSlider{Float64}(0.,10.,5.)
slider = Border(sl,"Slider")

imgout = ImageWidget()
imgw = Border(imgout, "Plot")

repl = Base.REPL.LineEditREPL(tty).mi
we = TerminalUI.WidgetREPL()
we.mi = repl
interface = Base.REPL.setup_interface(we,repl)

we.widget = w = TerminalUI.WidgetForREPL(interface,we)
run_repl(we)
b = WidgetStack([ slider, imgw, #RowLayout(gauges,[3 for i = 1:4],zeros(Int64,4),12),
    Border(w,"REPL")])
    #=Border(historyview,"History"),
    Border(,"REPL")])=#
import TerminalUI: children
TerminalUI.initialize!(b)
#TerminalUI.focus(children(children(w)[2])[2])
TerminalUI.focus(sl)
s = TerminalUI.DoubleBufferedTerminalScreen(Base.size(tty))

function resized(size)
    # Throw away all state, since we have no idea what the resize did
    #global s = TerminalUI.DoubleBufferedTerminalScreen(size)
    TerminalUI.redraw(s,b)
    #@show takebuf_string(TerminalUI.render(s))
    TerminalUI.render(STDOUT,s)
end

#lift(resized,TerminalUI.monitor_resize(tty))

function redraw(args...)
    TerminalUI.redraw(s,b)
    #@show takebuf_string(TerminalUI.render(s))
    TerminalUI.render(STDOUT,s)
end

function embed_here()
    try
        Base.reseteof(tty)
        Base.Terminals.raw!(tty,true)
        write(STDOUT,CSI,"?25l")
        write(STDOUT,CSI,"?6n")
        t = create_input_loop()
        wait(TerminalUI.curspos_condition)
        pos = TerminalUI.curspos
        global s
        s = TerminalUI.DoubleBufferedTerminalScreen((5,50); offset = pos, fullsize = Base.size(tty))
        redraw()
        wait(t)
        TerminalUI.afterembed(STDOUT,s)
    catch e
        rethrow(e)
    finally
        Base.Terminals.raw!(tty,false)
        write(STDOUT,CSI,"?25h")
    end
end

function pixel_for_cell_size(height, width)
    (height*34,width*14)
end

import VT100: Cell

function plot_to_imgout(p)
    buf = IOBuffer()
    size = (TerminalUI.optheight(imgout),TerminalUI.width(s))
    psize = pixel_for_cell_size(size...)
    # Macbook Pro
    dpi = 220.53
    png = PNG(buf,psize[2]/dpi*inch,psize[1]/dpi*inch; dpi = round(dpi))
    png.emit_on_finish = false
    #p = deepcopy(p)
    if isa(p,Gadfly.Plot)
        p.theme = Gadfly.Theme(
        panel_fill=color("black"), default_color=color("orange"), major_label_color=color("white"), minor_label_color=color("white"), key_label_color=color("white"), key_title_color=color("white")
       )
    end
    Compose.draw(png,p)
    imgout.data = takebuf_array(buf)
    if imgout.screenchar.content != 0
        push!(s.dropped, Int(imgout.screenchar.content) - 0xF0000)
    end
    imgout.screenchar = Cell(Char(0))
end

using Gadfly
lift(sl.val) do t
    plot_to_imgout(plot(x->sin(x+t),0,4,
        Scale.x_continuous(minvalue=0,maxvalue=4),
        Scale.y_continuous(minvalue=-1,maxvalue=1),
        Guide.xlabel(nothing),
        Guide.ylabel(nothing)))
    redraw()
end

#plot_to_imgout(compose(context(),rectangle(),fill("red")))

if !Base.isinteractive()
    try
        Base.reseteof(tty)
        Base.Terminals.raw!(tty,true)
        write(STDOUT,CSI,"?25l")
        write(STDOUT,CSI,"?1002h")
        write(STDOUT,CSI,"?1005h")
        lift(resized,TerminalUI.monitor_resize(tty))
        lift(redraw,TerminalUI.invalidated)
        #=lift(every(1.0)) do args...
            redraw()
            return true
        end=#
        wait(create_input_loop())
    finally
        Base.Terminals.raw!(tty,false)
        write(STDOUT,CSI,"?25h")
        write(STDOUT,CSI,"?1002l")
    end
end


text:
import Base: ==, *, /, +, -

"""
    LazyTensor(b1[, b2], indices, operators[, factor=1])

Lazy implementation of a tensor product of operators.

The suboperators are stored in the `operators` field. The `indices` field
specifies in which subsystem the corresponding operator lives. Note that these
must be sorted.
Additionally, a factor is stored in the `factor` field which allows for fast
multiplication with numbers.
"""
mutable struct LazyTensor{BL,BR,F,I,T} <: AbstractOperator{BL,BR}
    basis_l::BL
    basis_r::BR
    factor::F
    indices::I
    operators::T
    function LazyTensor(bl::BL, br::BR, indices::I, ops::T, factor::F=_default_factor(ops)) where {BL<:CompositeBasis,BR<:CompositeBasis,F,I,T<:Tuple}
        N = length(bl.bases)
        @assert N==length(br.bases)
        check_indices(N, indices)
        @assert length(indices) == length(ops)
        @assert issorted(indices)
        for n=1:length(indices)
            @assert isa(ops[n], AbstractOperator)
            @assert ops[n].basis_l == bl.bases[indices[n]]
            @assert ops[n].basis_r == br.bases[indices[n]]
        end
        F_ = promote_type(F,map(eltype, ops)...)
        factor_ = convert(F_, factor)
        new{BL,BR,F_,I,T}(bl, br, factor_, indices, ops)
    end
end
function LazyTensor(bl::CompositeBasis, br::CompositeBasis, indices, ops::Vector, factor=_default_factor(ops))
    Base.depwarn("LazyTensor(bl, br, indices, ops::Vector, factor) is deprecated, use LazyTensor(bl, br, indices, Tuple(ops), factor) instead.",
                :LazyTensor; force=true)
    return LazyTensor(bl,br,indices,Tuple(ops),factor)
end

function LazyTensor(basis_l::CompositeBasis, basis_r::Basis, indices::I, ops, factor::F=_default_factor(ops)) where {F,I}
    br = CompositeBasis(basis_r.shape, [basis_r])
    return LazyTensor(basis_l, br, indices, ops, factor)
end
function LazyTensor(basis_l::Basis, basis_r::CompositeBasis, indices::I, ops, factor::F=_default_factor(ops)) where {F,I}
    bl = CompositeBasis(basis_l.shape, [basis_l])
    return LazyTensor(bl, basis_r, indices, ops, factor)
end
function LazyTensor(basis_l::Basis, basis_r::Basis, indices::I, ops, factor::F=_default_factor(ops)) where {F,I}
    bl = CompositeBasis(basis_l.shape, [basis_l])
    br = CompositeBasis(basis_r.shape, [basis_r])
    return LazyTensor(bl, br, indices, ops, factor)
end
LazyTensor(op::T, factor) where {T<:LazyTensor} = LazyTensor(op.basis_l, op.basis_r, op.indices, op.operators, factor)
LazyTensor(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Integer, operator::T, factor=one(eltype(operator))) where T<:AbstractOperator = LazyTensor(basis_l, basis_r, [index], (operator,), factor)
LazyTensor(basis::Basis, index, operators, factor=_default_factor(operators)) = LazyTensor(basis, basis, index, operators, factor)

Base.copy(x::LazyTensor) = LazyTensor(x.basis_l, x.basis_r, copy(x.indices), Tuple(copy(op) for op in x.operators), x.factor)
Base.eltype(x::LazyTensor) = promote_type(eltype(x.factor), eltype.(x.operators)...)

function _default_factor(ops)
    Ts = map(eltype, ops)
    T = promote_type(Ts...)
    return one(T)
end
function _default_factor(op::T) where T<:AbstractOperator
    return one(eltype(op))
end

"""
    suboperator(op::LazyTensor, index)

Return the suboperator corresponding to the subsystem specified by `index`. Fails
if there is no corresponding operator (i.e. it would be an identity operater).
"""
suboperator(op::LazyTensor, index::Integer) = op.operators[findfirst(isequal(index), op.indices)]

"""
    suboperators(op::LazyTensor, index)

Return the suboperators corresponding to the subsystems specified by `indices`. Fails
if there is no corresponding operator (i.e. it would be an identity operater).
"""
suboperators(op::LazyTensor, indices) = [op.operators[[findfirst(isequal(i), op.indices) for i in indices]]...]

DenseOperator(op::LazyTensor) = op.factor*embed(op.basis_l, op.basis_r, op.indices, DenseOpType[DenseOperator(x) for x in op.operators])
SparseArrays.sparse(op::LazyTensor) = op.factor*embed(op.basis_l, op.basis_r, op.indices, SparseOpType[SparseOperator(x) for x in op.operators])

==(x::LazyTensor, y::LazyTensor) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.operators==y.operators && x.factor==y.factor

# Arithmetic operations
-(a::LazyTensor) = LazyTensor(a, -a.factor)

function +(a::LazyTensor{B1,B2}, b::LazyTensor{B1,B2}) where {B1,B2}
    if length(a.indices) == 1 && a.indices == b.indices
        op = a.operators[1] * a.factor + b.operators[1] * b.factor
        return LazyTensor(a.basis_l, a.basis_r, a.indices, (op,))
    end
    throw(ArgumentError("Addition of LazyTensor operators is only defined in case both operators act nontrivially on the same, single tensor factor."))
end

function -(a::LazyTensor{B1,B2}, b::LazyTensor{B1,B2}) where {B1,B2}
    if length(a.indices) == 1 && a.indices == b.indices
        op = a.operators[1] * a.factor - b.operators[1] * b.factor
        return LazyTensor(a.basis_l, a.basis_r, a.indices, (op,))
    end
    throw(ArgumentError("Subtraction of LazyTensor operators is only defined in case both operators act nontrivially on the same, single tensor factor."))
end

function *(a::LazyTensor{B1,B2}, b::LazyTensor{B2,B3}) where {B1,B2,B3}
    indices = sort(union(a.indices, b.indices))
    # ops = Vector{AbstractOperator}(undef, length(indices))
    function _prod(n)
        i = indices[n]
        in_a = i in a.indices
        in_b = i in b.indices
        if in_a && in_b
            return suboperator(a, i)*suboperator(b, i)
        elseif in_a
            a_i = suboperator(a, i)
            return a_i*identityoperator(typeof(a_i), b.basis_l.bases[i], b.basis_r.bases[i])
        elseif in_b
            b_i = suboperator(b, i)
            return identityoperator(typeof(b_i), a.basis_l.bases[i], a.basis_r.bases[i])*b_i
        end
    end
    ops = Tuple(_prod(n) for n ∈ 1:length(indices))
    return LazyTensor(a.basis_l, b.basis_r, indices, ops, a.factor*b.factor)
end
*(a::LazyTensor, b::Number) = LazyTensor(a, a.factor*b)
*(a::Number, b::LazyTensor) = LazyTensor(b, a*b.factor)
function *(a::LazyTensor{B1,B2}, b::DenseOpType{B2,B3}) where {B1,B2,B3}
    result = DenseOperator(a.basis_l, b.basis_r)
    mul!(result,a,b,complex(1.),complex(1.))
    result
end
function *(a::DenseOpType{B1,B2}, b::LazyTensor{B2,B3}) where {B1,B2,B3}
    result = DenseOperator(a.basis_l, b.basis_r)
    mul!(result,a,b,complex(1.),complex(1.))
    result
end

/(a::LazyTensor, b::Number) = LazyTensor(a, a.factor/b)

dagger(op::LazyTensor) = LazyTensor(op.basis_r, op.basis_l, op.indices, Tuple(dagger(x) for x in op.operators), conj(op.factor))

tensor(a::LazyTensor, b::LazyTensor) = LazyTensor(a.basis_l ⊗ b.basis_l, a.basis_r ⊗ b.basis_r, [a.indices; b.indices .+ length(a.basis_l.bases)], (a.operators..., b.operators...), a.factor*b.factor)

function tr(op::LazyTensor)
    b = basis(op)
    result = op.factor
    for i in 1:length(b.bases)
        if i in op.indices
            result *= tr(suboperator(op, i))
        else
            result *= length(b.bases[i])
        end
    end
    result
end

function ptrace(op::LazyTensor, indices)
    check_ptrace_arguments(op, indices)
    N = length(op.basis_l.shape)
    rank = N - length(indices)
    factor = op.factor
    for i in indices
        if i in op.indices
            factor *= tr(suboperator(op, i))
        else
            factor *= length(op.basis_l.bases[i])
        end
    end
    remaining_indices = remove(op.indices, indices)
    if rank==1 && length(remaining_indices)==1
        return factor * suboperator(op, remaining_indices[1])
    end
    b_l = ptrace(op.basis_l, indices)
    b_r = ptrace(op.basis_r, indices)
    if rank==1
        return factor * identityoperator(b_l, b_r)
    end
    ops = Tuple(suboperator(op, idx) for idx ∈ remaining_indices)
    LazyTensor(b_l, b_r, shiftremove(op.indices, indices), ops, factor)
end

normalize!(op::LazyTensor) = (op.factor /= tr(op); op)

function permutesystems(op::LazyTensor, perm)
    b_l = permutesystems(op.basis_l, perm)
    b_r = permutesystems(op.basis_r, perm)
    indices = [findfirst(isequal(i), perm) for i in op.indices]
    perm_ = sortperm(indices)
    LazyTensor(b_l, b_r, indices[perm_], op.operators[perm_], op.factor)
end

identityoperator(::Type{LazyTensor}, ::Type{T}, b1::Basis, b2::Basis) where T<:Number = LazyTensor(b1, b2, Int[], Tuple{}(), one(T))

# Recursively calculate result_{IK} = \\sum_J op_{IJ} h_{JK}
function _gemm_recursive_dense_lazy(i_k, N_k, K, J, val,
                        shape, strides_k, strides_j,
                        indices, h::LazyTensor,
                        op::Matrix, result::Matrix)
    if i_k > N_k
        for I=1:size(op, 1)
            result[I, K] += val*op[I, J]
        end
        return nothing
    end
    if i_k in indices
        h_i = suboperator(h, i_k)
        if isa(h_i, SparseOpPureType)
            h_i_data = h_i.data
            @inbounds for k=1:h_i_data.n
                K_ = K + strides_k[i_k]*(k-1)
                @inbounds for jptr=h_i_data.colptr[k]:h_i_data.colptr[k+1]-1
                    j = h_i_data.rowval[jptr]
                    J_ = J + strides_j[i_k]*(j-1)
                    val_ = val*h_i_data.nzval[jptr]
                    _gemm_recursive_dense_lazy(i_k+1, N_k, K_, J_, val_, shape, strides_k, strides_j, indices, h, op, result)
                end
            end
        elseif isa(h_i, DataOperator)
            h_i_data = h_i.data
            Nk = size(h_i_data, 2)
            Nj = size(h_i_data, 1)
            @inbounds for k=1:Nk
                K_ = K + strides_k[i_k]*(k-1)
                @inbounds for j=1:Nj
                    J_ = J + strides_j[i_k]*(j-1)
                    val_ = val*h_i_data[j,k]
                    _gemm_recursive_dense_lazy(i_k+1, N_k, K_, J_, val_, shape, strides_k, strides_j, indices, h, op, result)
                end
            end
        else
            throw(ArgumentError("gemm! of LazyTensor is not implemented for $(typeof(h_i))"))
        end
    else
        @inbounds for k=1:shape[i_k]
            K_ = K + strides_k[i_k]*(k-1)
            J_ = J + strides_j[i_k]*(k-1)
            _gemm_recursive_dense_lazy(i_k + 1, N_k, K_, J_, val, shape, strides_k, strides_j, indices, h, op, result)
        end
    end
end

# Recursively calculate result_{JI} = \\sum_K h_{JK} op_{KI}
function _gemm_recursive_lazy_dense(i_k, N_k, K, J, val,
                        shape, strides_k, strides_j,
                        indices, h::LazyTensor,
                        op::Matrix, result::Matrix)
    if i_k > N_k
        for I=1:size(op, 2)
            result[J, I] += val*op[K, I]
        end
        return nothing
    end
    if i_k in indices
        h_i = suboperator(h, i_k)
        if isa(h_i, SparseOpPureType) # TODO: adjoint sparse matrices
            h_i_data = h_i.data
            @inbounds for k=1:h_i_data.n
                K_ = K + strides_k[i_k]*(k-1)
                @inbounds for jptr=h_i_data.colptr[k]:h_i_data.colptr[k+1]-1
                    j = h_i_data.rowval[jptr]
                    J_ = J + strides_j[i_k]*(j-1)
                    val_ = val*h_i_data.nzval[jptr]
                    _gemm_recursive_lazy_dense(i_k+1, N_k, K_, J_, val_, shape, strides_k, strides_j, indices, h, op, result)
                end
            end
        elseif isa(h_i, DataOperator)
            h_i_data = h_i.data
            Nk = size(h_i_data, 2)
            Nj = size(h_i_data, 1)
            @inbounds for k=1:Nk
                K_ = K + strides_k[i_k]*(k-1)
                @inbounds for j=1:Nj
                    J_ = J + strides_j[i_k]*(j-1)
                    val_ = val*h_i_data[j,k]
                    _gemm_recursive_lazy_dense(i_k+1, N_k, K_, J_, val_, shape, strides_k, strides_j, indices, h, op, result)
                end
            end
        else
            throw(ArgumentError("gemm! of LazyTensor is not implemented for $(typeof(h_i))"))
        end
    else
        @inbounds for k=1:shape[i_k]
            K_ = K + strides_k[i_k]*(k-1)
            J_ = J + strides_j[i_k]*(k-1)
            _gemm_recursive_lazy_dense(i_k + 1, N_k, K_, J_, val, shape, strides_k, strides_j, indices, h, op, result)
        end
    end
end

function gemm(alpha, op::Matrix, h::LazyTensor, beta, result::Matrix)
    if iszero(beta)
        fill!(result, beta)
    elseif !isone(beta)
        rmul!(result, beta)
    end
    N_k = length(h.basis_r.bases)
    shape = [min(h.basis_l.shape[i], h.basis_r.shape[i]) for i=1:length(h.basis_l.shape)]
    strides_j = _strides(h.basis_l.shape)
    strides_k = _strides(h.basis_r.shape)
    _gemm_recursive_dense_lazy(1, N_k, 1, 1, alpha*h.factor, shape, strides_k, strides_j, h.indices, h, op, result)
end

function gemm(alpha, h::LazyTensor, op::Matrix, beta, result::Matrix)
    if iszero(beta)
        fill!(result, beta)
    elseif !isone(beta)
        rmul!(result, beta)
    end
    N_k = length(h.basis_l.bases)
    shape = [min(h.basis_l.shape[i], h.basis_r.shape[i]) for i=1:length(h.basis_l.shape)]
    strides_j = _strides(h.basis_l.shape)
    strides_k = _strides(h.basis_r.shape)
    _gemm_recursive_lazy_dense(1, N_k, 1, 1, alpha*h.factor, shape, strides_k, strides_j, h.indices, h, op, result)
end

mul!(result::DenseOpType{B1,B3},h::LazyTensor{B1,B2},op::DenseOpType{B2,B3},alpha,beta) where {B1,B2,B3} = (gemm(alpha, h, op.data, beta, result.data); result)
mul!(result::DenseOpType{B1,B3},op::DenseOpType{B1,B2},h::LazyTensor{B2,B3},alpha,beta) where {B1,B2,B3} = (gemm(alpha, op.data, h, beta, result.data); result)

function mul!(result::Ket{B1},a::LazyTensor{B1,B2},b::Ket{B2},alpha,beta) where {B1,B2}
    b_data = reshape(b.data, length(b.data), 1)
    result_data = reshape(result.data, length(result.data), 1)
    gemm(alpha, a, b_data, beta, result_data)
    result
end

function mul!(result::Bra{B2},a::Bra{B1},b::LazyTensor{B1,B2},alpha,beta) where {B1,B2}
    a_data = reshape(a.data, 1, length(a.data))
    result_data = reshape(result.data, 1, length(result.data))
    gemm(alpha, a_data, b, beta, result_data)
    result
end

function _tp_matmul_first!(result, a::AbstractMatrix, b, α::Number, β::Number)
    br = reshape(b, size(b, 1), :)
    result_r = reshape(result, size(a, 1), size(br, 2))
    mul!(result_r, a, br, α, β)
    result
end

function _tp_matmul_last!(result, a::AbstractMatrix, b, α::Number, β::Number)
    br = reshape(b, :, size(b, ndims(b)))
    result_r = reshape(result, (size(br, 1), size(a, 1)))
    mul!(result_r, br, transpose(a), α, β)
    result
end

function _tp_matmul_mid!(result, a::AbstractMatrix, loc::Integer, b, α::Number, β::Number)
    shp_b_1 = 1
    for i in 1:loc-1
        shp_b_1 *= size(b,i)
    end
    shp_b_3 = 1
    for i in loc+1:ndims(b)
        shp_b_3 *= size(b,i)
    end

    # TODO: Perhaps we should avoid reshaping here... should be possible to infer
    # contraction index tuple sizes
    br = reshape(b, shp_b_1, size(b, loc), shp_b_3)
    result_r = reshape(result, shp_b_1, size(a, 1), shp_b_3)

    # If using BLAS, this will generally require intermediate storage
    # TensorOperations will use its own cache for these.

    # In principle, we can use the @tensor macro here and have TensorOperations
    # take care of the cache symbols. However, this seems to be broken...?
    #@tensor result_r[a,b,c] = α * a[b,x] * br[a,x,c] + β * result_r[a,b,c]

    # these are used to identify objects in the cache (if used)
    syms = (:_tp_matmul_a, :_tp_matmul_b, :_tp_matmul_c)

    oindA, cindA, oindB, cindB, indCinoAB = TensorOperations.contract_indices(
      (-2, 1), (-1, 1, -3), (-1, -2, -3))

    TensorOperations.contract!(
      α, a, :N, br, :N, β, result_r, oindA, cindA, oindB, cindB, indCinoAB, syms)

    result
end

function _tp_matmul!(result, a::AbstractMatrix, loc::Integer, b, α::Number, β::Number)
    """Apply a matrix `a` to one tensor factor of a tensor `b`.

    If β is nonzero, add to β times `result`. In other words, we do:

    result = α * a * b + β * result

    Parameters:
        result: Array to hold the output tensor.
        a: Matrix to apply.
        loc: Index of the dimension of `b` to which `a` should be applied.
        b: Array containing the tensor to which `a` should be applied.
        α: Factor multiplying `a`.
        β: Factor multiplying `result`.
    Returns:
        The modified `result`.
    """
    # Use GEMM directly where possible, otherwise TensorOperations
    if loc == 1
        return _tp_matmul_first!(result, a, b, α, β)
    elseif loc == ndims(b)
        return _tp_matmul_last!(result, a, b, α, β)
    end
    _tp_matmul_mid!(result, a, loc, b, α, β)
end

function _tp_sum_get_tmp(op::AbstractMatrix{T}, loc::Integer, arr::AbstractArray{T,N}, sym) where {T,N}
    shp = ntuple(i -> i == loc ? size(op,1) : size(arr,i), N)
    if TensorOperations.use_cache()
        tmp::Array{T,N} = get!(TensorOperations.cache, (sym, taskid(), T, shp)) do
            Array{T,N}(undef, shp...)
        end
    else
        tmp = Array{T,N}(undef, shp...)
    end

    tmp
end

#Apply a tensor product of operators to a vector.
function _tp_sum_matmul!(result_data, tp_ops, iso_ops, b_data, alpha, beta)
    if iso_ops === nothing
        n_ops = length(tp_ops[1])
        ops = zip(tp_ops...)
    else
        n_ops = length(tp_ops[1]) + length(iso_ops[1])
        ops = Iterators.flatten((zip(tp_ops...), zip(iso_ops...)))
    end

    # TODO: Perhaps replace with a single for loop and branch inside?
    if n_ops == 0
        result_data .= alpha .* b_data + beta .* result_data
    elseif n_ops == 1
        # Can add directly to the output array.
        _tp_matmul!(result_data, first(ops)..., b_data, alpha, beta)
    elseif n_ops == 2
        # One temporary vector needed.
        op1, istate = iterate(ops)
        tmp = _tp_sum_get_tmp(op1..., b_data, :_tp_sum_matmul_tmp1)
        _tp_matmul!(tmp, op1..., b_data, alpha, zero(beta))

        op2, istate = iterate(ops, istate)
        _tp_matmul!(result_data, op2..., tmp, one(alpha), beta)
    else
        # At least two temporary vectors needed.
        # Symbols identify computation stages in the TensorOperations cache.
        sym1 = :_tp_sum_matmul_tmp1
        sym2 = :_tp_sum_matmul_tmp2

        op1, istate = iterate(ops)
        tmp1 = _tp_sum_get_tmp(op1..., b_data, sym1)
        _tp_matmul!(tmp1, op1..., b_data, alpha, zero(beta))

        next = iterate(ops, istate)
        for _ in 2:n_ops-1
            op, istate = next
            tmp2 = _tp_sum_get_tmp(op..., tmp1, sym2)
            _tp_matmul!(tmp2, op..., tmp1, one(alpha), zero(beta))
            tmp1, tmp2 = tmp2, tmp1
            sym1, sym2 = sym2, sym1
            next = iterate(ops, istate)
        end

        op, istate = next
        _tp_matmul!(result_data, op..., tmp1, one(alpha), beta)
    end

    result_data
end

# Represents a rectangular "identity" matrix of ones along the diagonal.
struct _SimpleIsometry{I<:Integer}  # We have no need to subtype AbstractMatrix
    shape::Tuple{I,I}
    function _SimpleIsometry(d1::I, d2::I) where {I<:Integer}
        new{I}((d1, d2))
    end
end
Base.size(A::_SimpleIsometry) = A.shape
Base.size(A::_SimpleIsometry, i) = A.shape[i]

function _tp_matmul!(result, a::_SimpleIsometry, loc::Integer, b, α::Number, β::Number)
    shp_b_1 = 1
    for i in 1:loc-1
        shp_b_1 *= size(b,i)
    end
    shp_b_3 = 1
    for i in loc+1:ndims(b)
        shp_b_3 *= size(b,i)
    end

    @assert size(b, loc) == size(a, 2)

    br = Base.ReshapedArray(b, (shp_b_1, size(b, loc), shp_b_3), ())
    result_r = Base.ReshapedArray(result, (shp_b_1, size(a, 1), shp_b_3), ())

    d = min(size(a)...)

    if β != 0
        rmul!(result, β)
        result_r[:, 1:d, :] .+= α .* br[:, 1:d, :]
    else
        result_r[:, 1:d, :] .= α .* br[:, 1:d, :]
        result_r[:, d+1:end, :] .= zero(eltype(result))
    end

    result
end

function _explicit_isometries(used_indices, bl::Basis, br::Basis, shift=0)
    indices = Set(used_indices)
    isos = nothing
    iso_inds = nothing
    for (i, (sl, sr)) in enumerate(zip(_comp_size(bl), _comp_size(br)))
        if (sl != sr) && !(i + shift in indices)
            if isos === nothing
                isos = [_SimpleIsometry(sl, sr)]
                iso_inds = [i + shift]
            else
                push!(isos, _SimpleIsometry(sl, sr))
                push!(iso_inds, i + shift)
            end
        end
    end
    if isos === nothing
        return nothing
    end
    isos, iso_inds
end

# To get the shape of a CompositeBasis with number of dims inferrable at compile-time
_comp_size(b::CompositeBasis) = tuple((length(b_) for b_ in b.bases)...)
_comp_size(b::Basis) = (length(b),)

function mul!(result::Ket{B1}, a::LazyTensor{B1,B2,F,I,T}, b::Ket{B2}, alpha, beta) where {B1,B2, F,I,T<:Tuple{Vararg{DenseOpType}}}
    # We reshape here so that we have the proper shape information for the
    # tensor contraction later on. Using ReshapedArray vs. reshape() avoids
    # an allocation.
    b_data = Base.ReshapedArray(b.data, _comp_size(basis(b)), ())
    result_data = Base.ReshapedArray(result.data, _comp_size(basis(result)), ())

    tp_ops = (tuple((op.data for op in a.operators)...), a.indices)
    iso_ops = _explicit_isometries(a.indices, a.basis_l, a.basis_r)
    _tp_sum_matmul!(result_data, tp_ops, iso_ops, b_data, alpha * a.factor, beta)

    result
end

function mul!(result::Bra{B2}, a::Bra{B1}, b::LazyTensor{B1,B2,F,I,T}, alpha, beta) where {B1,B2, F,I,T<:Tuple{Vararg{DenseOpType}}}
    a_data = Base.ReshapedArray(a.data, _comp_size(basis(a)), ())
    result_data = Base.ReshapedArray(result.data, _comp_size(basis(result)), ())

    tp_ops = (tuple((transpose(op.data) for op in b.operators)...), b.indices)
    iso_ops = _explicit_isometries(b.indices, b.basis_r, b.basis_l)
    _tp_sum_matmul!(result_data, tp_ops, iso_ops, a_data, alpha * b.factor, beta)

    result
end

function mul!(result::DenseOpType{B1,B2}, a::LazyTensor{B1,B1,F,I,T}, b::DenseOpType{B1,B2}, alpha, beta) where {B1,B2, F,I,T<:Tuple{Vararg{DenseOpType}}}
    b_data = Base.ReshapedArray(b.data, (_comp_size(b.basis_l)..., _comp_size(b.basis_r)...), ())
    result_data = Base.ReshapedArray(result.data, (_comp_size(result.basis_l)..., _comp_size(result.basis_r)...), ())

    tp_ops = (tuple((op.data for op in a.operators)...), a.indices)
    iso_ops = _explicit_isometries(a.indices, a.basis_l, a.basis_r)
    _tp_sum_matmul!(result_data, tp_ops, iso_ops, b_data, alpha * a.factor, beta)

    result
end

function mul!(result::DenseOpType{B1,B3}, a::DenseOpType{B1,B2}, b::LazyTensor{B2,B3,F,I,T}, alpha, beta) where {B1,B2,B3, F,I,T<:Tuple{Vararg{DenseOpType}}}
    a_data = Base.ReshapedArray(a.data, (_comp_size(a.basis_l)..., _comp_size(a.basis_r)...), ())
    result_data = Base.ReshapedArray(result.data, (_comp_size(result.basis_l)..., _comp_size(result.basis_r)...), ())

    shft = length(a.basis_l.shape)  # b must be applied to the "B2" side of a
    tp_ops = (tuple((transpose(op.data) for op in b.operators)...), tuple((i + shft for i in b.indices)...))
    iso_ops = _explicit_isometries(tp_ops[2], b.basis_r, b.basis_l, shft)
    _tp_sum_matmul!(result_data, tp_ops, iso_ops, a_data, alpha * b.factor, beta)

    result
end


text:
export Gradient

################################################################################
######   Preconditioning algorithm definition (structure holding params)  ######
################################################################################
"""
    Gradient()

Algorithm for descending along the steepest gradient with SGD-based optimizers.
"""
struct Gradient <: Algorithm end

################################################################################
######  Structure holding information computed at the end of a sampling   ######
################################################################################
"""
    GradientEvaluation <: EvaluatedAlgorithm

The GradientEvaluation structure holds the evaluation of the Network, L, the
vector of generalized forces F acting on it.
"""
mutable struct GradientEvaluation{TL,TF} <: EvaluatedAlgorithm
    L::TL
    F::TF

    # Individual values to compute statistical correlators
    LVals::Vector
end

function GradientEvaluation(net::NeuralNetwork)
    wt = grad_cache(net)
    WT = weight_type(net)
    T = out_type(net)

    F = Tuple([zeros(WT,size(w)) for w=wt.tuple_all_weights])

    GradientEvaluation(zero(T),
                       F,
                       Vector{T}())
end

EvaluatedNetwork(alg::Gradient, net) =
    GradientEvaluation(weights(net))

# Utility method utilised to accumulate results on a single variable
function add!(acc::GradientEvaluation, o::GradientEvaluation)
   acc.L  += o.L
   acc.F .+= o.F

   append!(acc.LVals, o.LVals)
end

function precondition!(∇x, params::Gradient, data::GradientEvaluation, args...)
    for (Δw, F) = zip(∇x, data.F)
        Δw .= F
    end
    true
end


text:
## a list of reserved names from PyCall.jl
const reserved = Set(["while", "if", "for", "try", "return", "break",
    "continue", "function", "macro", "quote", "let", "local",
    "global", "const", "abstract", "typealias", "type",
    "bitstype", "immutable", "ccall", "do", "module",
    "baremodule", "using", "import", "export", "importall",
    "false", "true"])

function rwrap(pkg::Compat.String, s::Symbol)
    reval("library($pkg)")
    members = rcopy("ls('package:$pkg')")
    filter!(x -> !(x in reserved), members)
    m = Module(s)
    consts = [:(const $(Symbol(x)) = rcall(Symbol("::"),$(QuoteNode(Symbol(pkg))),$(QuoteNode(Symbol(x))))) for x in members]
    id = Expr(:(=), :__package__, pkg)
    exports = [Symbol(x) for x in members]
    s in exports && error("$pkg has a function with the same name as $(pkg), use `@rimport $pkg as ...` instead.")
    eval(m, Expr(:toplevel, :(import RCall.rcall), consts..., Expr(:export, exports...), id, :(__exports__ = $exports)))
    m
end

"Import an R Package as a Julia module. You can also use classic Python syntax to make an alias: `@rimport *module-name* as *shorthand*`"
macro rimport(x, args...)
    if length(args)==2 && args[1] == :as
        m = args[2]
    elseif length(args)==0
        m = x
    else
        throw(ArgumentError("invalid import syntax."))
    end
    pkg = string(x)
    sym = Expr(:quote, m)
    quote
        if !isdefined($sym)
            const $(esc(m)) = rwrap($pkg, $sym)
            nothing
        elseif typeof($(esc(m))) <: Module &&
                    :__package__ in names($(esc(m)), true) &&
                    $(esc(m)).__package__ == $pkg
            nothing
        else
            error("$($sym) already exists!")
            nothing
        end
    end
end

"""
Load all exported functions/objects of a R package to the current module.
"""
macro rlibrary(x)
    pkg = Expr(:quote, x)
    quote
        reval("library($($pkg))")
        members = rcopy("ls('package:$($pkg)')")
        filter!(x -> !(x in reserved), members)
        for m in members
            sym = Symbol(m)
            eval(current_module(), Expr(
                    :(=),
                    sym,
                    Expr(:call, :rcall, QuoteNode(Symbol("::")), QuoteNode($pkg), QuoteNode(sym))
                )
            )
        end
    end
end


text:
mutable struct Advance{T} <: State{T}
    value::T
    t::T
    Δt::T
end

Advance(; init=nothing, step=nothing, unit, _type, _...) = begin
    U = value(unit)
    T = valuetype(_type, U)
    t = isnothing(init) ? zero(T) : unitfy(value(init), U)
    Δt = isnothing(step) ? oneunit(T) : unitfy(value(step), U)
    #T = promote_type(typeof(t), typeof(Δt))
    Advance{T}(t, t, Δt)
end

constructortags(::Val{:Advance}) = (:init, :step, :unit)

genvartype(v::VarInfo, ::Val{:Advance}; V, _...) = @q Advance{$V}

gendefault(v::VarInfo, ::Val{:Advance}) = missing

genupdate(v::VarInfo, ::Val{:Advance}, ::MainStep; kw...) = begin
    @gensym s t
    @q let $s = $(symstate(v)),
           $t = $s.t
        $s.t += $s.Δt
        $C.store!($s, $t)
    end
end


text:
using Juno
import Juno: Tree, icon

media(SyncPlot, Media.Plot)
media(Plot, Media.Plot)

Media.render(::Juno.PlotPane, plot::SyncPlot) = display_blink(plot)

@render Juno.PlotPane plot::Plot SyncPlot(plot)

@render Juno.Editor p::SyncPlot p.plot

@render Juno.Editor p::Plot begin
  Tree(icon("graph"), [p.data, p.layout])
end


text:
include("RCF/kummer_extensions.jl")
include("RCF/class_fields.jl")
include("RCF/cyclo.jl")
include("RCF/rcf.jl")
include("RCF/autos.jl")
include("RCF/conductor.jl")
include("RCF/rcf_subfield.jl")


text:
module MCMCStorage

export Chains, StanCSV

include("chains.jl")
include("stan_csv.jl")

end # module


text:
using Falcons
using Healpix
using NPZ
using PyCall
using Base.Threads
using StaticArrays
using Formatting
using LinearAlgebra
using HDF5

function save_psiDB(save_dir, idx, db)
    h5open(save_dir * "idx=$idx.h5", "w") do file
        write(file,"map_div", db[1])
        write(file,"idx", db[2])
        for i in eachindex(db[3])
            #println(i)
            write(file,"/psi/$i", db[3][i])
        end
    end
end

hp = pyimport("healpy")
np = pyimport("numpy")
pd = pyimport("pandas")
glob = pyimport("glob")
re = pyimport("re")
os = pyimport("os")

day = 60 * 60 * 24
year = day * 365

ss = gen_ScanningStrategy()
ss.nside = 512
ss.duration = 3year #[sec]
ss.sampling_rate = 19 #[Hz]
ss.FP_theta = [0]#mft_focalplane_pointings[:,3]
ss.FP_phi = [0]#mft_focalplane_pointings[:,4].+90
ss.alpha = 45 #[degree]
ss.beta = 50 #[degree]
ss.prec_rpm = period2rpm(192.348) #[sec]
ss.spin_rpm = 0.05 #[rpm]
ss.hwp_rpm = 0 #[rpm]
ss.start_point = "equator"

idx = parse(Int64, ARGS[1])

dir = "/group/cmb/litebird/usr/ytakase/psi_db/"
det_pixel = "3_30/"#"$ss.FP_theta[1]_$ss.FP_phi[1]"#将来はピクセルに置き換えるべき
save_dir = dir * det_pixel

if os.path.exists(save_dir) == false
    os.makedirs(save_dir)
end

db = get_psiDataBase(ss, division=32*10, idx=idx, map_div=3);

@show sizeof(db)

save_psiDB(save_dir, idx, db)
#np.savez_compressed(save_dir*det_pixel*"idx=$idx", psi=psi_db)


text:
function MartinStrahlen(Dim,x0,x1,y0,y1,N,JJ,nx,ny,verschx,verschy)

laengenmatrix = zeros(1,nx*ny);
theta1 = 0.00001; #Referenzwinkel
Nullzeile = 1;
alpha = zeros(JJ,1);

for i=1:JJ
    theta= theta1+((i-1)/JJ)*pi;
    #omega[i] = theta;
    alpha[i] = theta*180/pi;
    X,Y,Deltas = createStartingPositionsmeinsneu(x0,x1,y0,y1,N,nx,ny,theta,verschx,verschy);
    vx = cos(theta);
    vy = sin(theta);
    laengenmatrixTeil, SchnittVektorsort,Nullzeile = Schnittpunkte(Dim,i,X,Y,x0,y0,x1,y1,N,JJ,nx,ny,vx,vy,theta,Nullzeile);
    if laengenmatrix == zeros(1,nx*ny)
        laengenmatrix = laengenmatrixTeil;
    else
        laengenmatrix = [laengenmatrix; laengenmatrixTeil];
    end
end
return laengenmatrix;
end


text:
export BasicDQNLearner

"""
    BasicDQNLearner(;kwargs...)

See paper: [Playing Atari with Deep Reinforcement Learning](https://www.cs.toronto.edu/~vmnih/docs/dqn.pdf)

This is the very basic implementation of DQN. Compared to the traditional Q learning, the only difference is that,
in the updating step it uses a batch of transitions sampled from an experience buffer instead of current transition.
And the `approximator` is usually a [`NeuralNetworkApproximator`](@ref).
You can start from this implementation to understand how everything is organized and how to write your own customized algorithm.

# Keywords

- `approximator`::[`AbstractApproximator`](@ref): used to get Q-values of a state.
- `loss_func`: the loss function to use.
- `γ::Float32=0.99f0`: discount rate.
- `batch_size::Int=32`
- `min_replay_history::Int=32`: number of transitions that should be experienced before updating the `approximator`.
- `rng=Random.GLOBAL_RNG`
"""
mutable struct BasicDQNLearner{Q,F,R} <: AbstractLearner
    approximator::Q
    loss_func::F
    γ::Float32
    sampler::BatchSampler
    min_replay_history::Int
    rng::R
    # for debugging
    loss::Float32
end

Flux.functor(x::BasicDQNLearner) = (Q = x.approximator,), y -> begin
    x = @set x.approximator = y.Q
    x
end

(learner::BasicDQNLearner)(env) =
    env |>
    get_state |>
    x -> send_to_device(device(learner), x) |>
    learner.approximator |>
    send_to_host

function BasicDQNLearner(;
    approximator::Q,
    loss_func::F = huber_loss,
    γ = 0.99f0,
    batch_size = 32,
    min_replay_history = 32,
    rng = Random.GLOBAL_RNG,
) where {Q,F}
    BasicDQNLearner{Q,F,typeof(rng)}(
        approximator,
        loss_func,
        γ,
        BatchSampler{SARTS}(batch_size),
        min_replay_history,
        rng,
        0.0,
    )
end

function RLBase.update!(learner::BasicDQNLearner, traj::AbstractTrajectory)
    if length(traj) >= learner.min_replay_history
        inds, batch = sample(learner.rng, traj, learner.sampler)
        update!(learner, batch)
    end
end

function RLBase.update!(learner::BasicDQNLearner, batch::NamedTuple{SARTS})

    Q = learner.approximator
    γ = learner.γ
    loss_func = learner.loss_func

    s, a, r, t, s′ = send_to_device(device(Q), batch)
    a = CartesianIndex.(a, 1:length(a))

    gs = gradient(params(Q)) do
        q = Q(s)[a]
        q′ = vec(maximum(Q(s′); dims = 1))
        G = r .+ γ .* (1 .- t) .* q′
        loss = loss_func(G, q)
        ignore() do
            learner.loss = loss
        end
        loss
    end

    update!(Q, gs)
end


text:
##file to include in ComplexModeler.jl which contains all files associated to GOC

include("GOCGenerator.jl")
include("GOCLoad.jl")
include("GOCShunt.jl")
include("GOCVolt.jl")
include("GOCLineπ_notransformer.jl")
include("GOCLineπ_withtransformer.jl")
include("GOCNullImpedance_notransformer.jl")
include("GOCNullImpedance_withtransformer.jl")
include("GOC_read.jl")
include("GOC_read_solutions.jl")
include("GOCcheck_feasibility.jl")
include("GOC_read_raw.jl")
include("GOC_write_solutions.jl")


text:
# Example for HINT layer (Kruse et al, 2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test

# Input
nx = 32
ny = 32
n_channel = 32
n_hidden = 64
batchsize = 2

# Input image
X = glorot_uniform(nx, ny, n_channel, batchsize)

###################################################################################################
# Basic HINT layer

# Create HINT layer w/o logdet
HL1 = CouplingLayerHINT(n_channel, n_hidden; logdet=false, permute="none")

# Call forward and inverse
Y = HL1.forward(X)
X_ = HL1.inverse(Y)

@test isapprox(norm(X_ - X)/norm(X), 0f0; atol=1f-5)

###################################################################################################
# Basic HINT layer w/ logdet

# HINT layer with logdet
HL2 = CouplingLayerHINT(n_channel, n_hidden; logdet=true, permute="none")

Y, logdet = HL2.forward(X)
X_ = HL2.inverse(Y)

@test isapprox(norm(X_ - X)/norm(X), 0f0; atol=1f-5)

###################################################################################################
# Basic HINT layer w/ logdet and permutation

# HINT layer with permutation: set to "none", "lower" or "full"
HL3 = CouplingLayerHINT(n_channel, n_hidden; logdet=true, permute="lower")

Y, logdet = HL3.forward(X)
X_ = HL3.inverse(Y)

@test isapprox(norm(X_ - X)/norm(X), 0f0; atol=1f-5)