text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

#=
function parse_commandline_error_handler(settings::ArgParseSettings, err, err_code::Int = 1)
	help_string =
		"execute one of the following to get the help instructions:\n"*
		"\t- from the CLI: \t\tjulia $(FILE_NAME) --help\n"*
		"\t- from the CLI: \t\t./$FILE_NAME --help\n"*
		"\t- inside the Julia REPL: \tinclude(\"$(FILE_NAME)\n)"

    if occursin("no command given", err.text)
     	println("\nERROR: no input command were given.\n"*help_string)
	elseif occursin("unknown command", err.text)
		println(
			"\nERROR, $(err.text)\n"*
			usage_string(settings)*"\n\n"*
			help_string)
    elseif occursin("out of range input for input_file:", err.text)
     	println(stderr, "input_file is not a PFM file")
	elseif  occursin("too many", err.text) &&  occursin("usage: $(FILE_NAME) demo", usage_string(settings))
		ArgParse.show_help(settings)
	else
     	println(stderr, err.text, "\n" , usage_string(settings))
		#ArgParse.show_help(settings)
    	end

	throw(ArgumentError("error code : $err_code"))
	#exit(err_code)
	#return err_code
end
=#

function ArgParse_command_line(arguments)
	s = ArgParseSettings()

	s.description = "Raytracer for the generation of photorealistic images in Julia."
	#s.exc_handler = parse_commandline_error_handler
	#s.version = @project_version

	@add_arg_table! s begin
		"demo"
			action = :command
			help = "create a standard image that checks the correct behaviour of the program"
		"demo_animation"
			action = :command
			help = "create an animation of a 360 degree rotation of the demo image"
		"tonemapping"
			action = :command
			help = "apply tone mapping to a pfm image and save it as a ldr file"
		"render"
			action = :command
			help = "render an image from a file"
		"animation"
			action = :command
			help = "render an animation from a file"
	end

	#### TONE MAPPING ###################################################################92

	s["tonemapping"].description = "Apply tone mapping to a pfm image and save it as a ldr file"

	add_arg_group!(s["tonemapping"], "tonemapping filenames");
	@add_arg_table! s["tonemapping"] begin
		"infile"
			help = "path to input file, it must be a PFM file"
			#range_tester = input -> (typeof(query(input))<:File{format"PFM"})
			required = true
		"outfile"
			help = "output file name"
			required = true
	end

	add_arg_group!(s["tonemapping"], "tonemapping settings");
	@add_arg_table! s["tonemapping"] begin
		"--normalization", "-a"
			help = "scaling factor for the normalization process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.18
		"--gamma", "-g"
			help = "gamma value for the tone mapping process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 1.27
          "--avg_lum"
			help = "average luminosity of the resulting pfm image; if =0 (dafault value)"*
                    "it will be calculated automatically with the built-in avg_lum function"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.0
		"--ONLY_FOR_TESTS"
			help = "only for testing flag; do not use it"
			action = :store_true
	end

	#### DEMO ###########################################################################92

	s["demo"].description =
		"""Creates a demo image with the specified options.\n"""*
		"""There are two possible demo image "world" to be """*
		"""rendered, specified through the input string `--world-type`.\n\n"""*
		"""The `type=="A"` demo image world consist in a set of 10 spheres of equal radius 0.1:"""*
		"""8 spheres are placed at the verteces of a cube of side 1.0, one in the center of"""*
		"""the lower cube face and the last one in the center of the left cube face.\n\n"""*
		"""The `type=="B"` demo image world consists in a checked x-y plane, a blue opaque"""*
		"""sphere, a red reflecting sphere, and a green oblique reflecting plane, all"""*
		"""inside a giant emetting sphere.\n\n"""*
		"""The creation of the demo image has the objective to check the correct behaviour of
		the rendering software, specifically the orientation upside-down and left-right."""

	add_arg_group!(s["demo"], "render options");
	@add_arg_table! s["demo"] begin
		"--camera_type"
			help = "option for the camera type:\n"*
	    				"ort -> Orthogonal camera, per -> Perspective camera"
          	arg_type = String
			default = "per"
			range_tester = input -> check_is_one_of(input, CAMERAS)
		"--camera_position"
          	help = "camera position in the scene as '[X,Y,Z]'"
          	arg_type = String
          	default = "[-1,0,0]"
          	range_tester = check_is_vector
		"--alpha"
			help = "angle of view around z-axis, in degrees"
			arg_type = Float64
			default = 0.
		"--width"
			help = "pixel number on the width of the resulting demo image."
			default = 640
			range_tester = check_is_even_uint64
		"--height"
			help = "pixel number on the height of the resulting demo image."
			default = 480
			range_tester = check_is_even_uint64
          "--normalization", "-a"
			help = "scaling factor for the normalization process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.18
		"--gamma", "-g"
			help = "gamma value for the tone mapping process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 1.27
          "--avg_lum"
			help = "average luminosity of the resulting pfm image; if =0 (dafault value)"*
                    "it will be calculated automatically with the built-in avg_lum function"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.0
     	"--samples_per_pixel"
			help = "Number of samples per pixel for the antialiasing algorithm\n"*
					"It must be an integer perfect square, i.e. 0,1,4,9,16,...\n"*
					"If =0 (default value), antialiasing does not occurs."
     		default = 0
			range_tester = check_is_square
		"--world_type"
          	help = "world type to be rendered; valid values: $(DEMO_WORLD_TYPES)"
          	arg_type = String
          	default = "A"
          	range_tester = input -> check_is_one_of(input, DEMO_WORLD_TYPES)
		"--set_pfm_name"
			help = "name of the pfm file to be saved"
			nargs = '?'
			arg_type = String
			default = "demo.pfm"
			constant = "demo.pfm"
		"--set_png_name"
			help = "name of the png file to be saved"
			nargs = '?'
			arg_type = String
			default = "demo.png"
			constant = "demo.png"
		"--ONLY_FOR_TESTS"
			help = "only for testing flag; do not use it"
			action = :store_true
	end

	add_arg_group!(s["demo"], "renderer to be used");
	@add_arg_table! s["demo"] begin
		"onoff"
			action = :command
			help = "onoff renderer"
			#dest_name = "onoff"
		"flat"
			action = :command
			help = "flat-renderer"
			#dest_name = "flat"
		"pathtracer"
			action = :command
			help = "path tracing renderer"
			#dest_name = "pathtracer"
		"pointlight"
			action = :command
			help = "point-light tracing renderer"
			#dest_name = "pointlight"
	end

	add_arg_group!(s["demo"]["onoff"], "onoff renderer options");
	@add_arg_table! s["demo"]["onoff"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--color"
			help = "hit color specified as '<R,G,B>' components. Example: --ambient_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
	end

	add_arg_group!(s["demo"]["flat"], "flat renderer options");
	@add_arg_table! s["demo"]["flat"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
	end

	add_arg_group!(s["demo"]["pathtracer"], "pathtracing renderer options");
	@add_arg_table! s["demo"]["pathtracer"] begin
		"--init_state"
    			help = "Initial seed for the random number generator (positive integer number)."
    			default = 45
			range_tester = check_is_uint64
    		"--init_seq"
    			help = "Identifier of the sequence produced by the "*
			    "random number generator (positive integer number)."
    			default = 54
			range_tester = check_is_uint64
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--num_of_rays"
			help = "number of `Ray`s generated for each integral evaluation"
			default = 10
			range_tester = check_is_uint64
		"--max_depth"
			help = "maximal number recursive integrations"
			default = 3
			range_tester = check_is_uint64
		"--russian_roulette_limit"
			help = "depth at whitch the Russian Roulette algorithm begins"
			default = 2
			range_tester = check_is_uint64
	end

	add_arg_group!(s["demo"]["pointlight"], "pointlight renderer options");
	@add_arg_table! s["demo"]["pointlight"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--ambient_color"
			help = "ambient color specified as '<R,G,B>' components. Example: --ambient_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--dark_parameter"
			help = "percentage of the retuned hit point color if it is not directly visible from a point-light source"
			default = 0.05
			range_tester = check_is_positive
	end

	#### DEMO_ANIMATION #################################################################92

	s["demo_animation"].description = "creates an animation of a 360 degree rotation around"*
								"vertical axis of the demo image."

	add_arg_group!(s["demo_animation"], "demo_animation settings");
	@add_arg_table! s["demo_animation"] begin
		"--camera_type"
			help = "option for the camera type:\n"*
	    				"ort -> Orthogonal camera, per -> Perspective camera"
          	arg_type = String
			default = "per"
			range_tester = input -> check_is_one_of(input, CAMERAS)
		"--camera_position"
          	help = "camera position in the scene as '[X,Y,Z]'"
          	arg_type = String
          	default = "[-1,0,0]"
          	range_tester = check_is_vector
		"--width"
			help = "pixel number on the width of the resulting demo image."
			default = 640
			range_tester = check_is_even_uint64
		"--height"
			help = "pixel number on the height of the resulting demo image."
			default = 480
			range_tester = check_is_even_uint64
          "--normalization", "-a"
			help = "scaling factor for the normalization process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.18
		"--gamma", "-g"
			help = "gamma value for the tone mapping process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 1.27
          "--avg_lum"
			help = "average luminosity of the resulting pfm image; if =0 (dafault value)"*
                    "it will be calculated automatically with the built-in avg_lum function"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.0
     	"--samples_per_pixel"
			help = "Number of samples per pixel for the antialiasing algorithm\n"*
					"It must be an integer perfect square, i.e. 0,1,4,9,16,...\n"*
					"If =0 (default value), antialiasing does not occurs."
     		default = 0
			range_tester = check_is_square
		"--world_type"
          	help = "world type to be rendered; valid values: $(DEMO_WORLD_TYPES)"
          	arg_type = String
          	default = "A"
          	range_tester = input -> check_is_one_of(input, DEMO_WORLD_TYPES)
		"--set_anim_name"
			help = "name of the animation file to be saved"
			nargs = '?'
			arg_type = String
			default = "demo_animation.mp4"
			constant = "demo_animation.mp4"
		"--ONLY_FOR_TESTS"
			help = "only for testing flag; do not use it"
			action = :store_true
	end

	add_arg_group!(s["demo_animation"], "renderer to be used");
	@add_arg_table! s["demo_animation"] begin
		"onoff"
			action = :command
			help = "onoff renderer"
			#dest_name = "onoff"
		"flat"
			action = :command
			help = "flat-renderer"
			#dest_name = "flat"
		"pathtracer"
			action = :command
			help = "path tracing renderer"
			#dest_name = "pathtracer"
		"pointlight"
			action = :command
			help = "point-light tracing renderer"
			#dest_name = "pointlight"
	end

	add_arg_group!(s["demo_animation"]["onoff"], "onoff renderer options");
	@add_arg_table! s["demo_animation"]["onoff"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--color"
			help = "hit color specified as '<R,G,B>' components. Example: --ambient_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
	end

	add_arg_group!(s["demo_animation"]["flat"], "flat renderer options");
	@add_arg_table! s["demo_animation"]["flat"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
	end

	add_arg_group!(s["demo_animation"]["pathtracer"], "pathtracing renderer options");
	@add_arg_table! s["demo_animation"]["pathtracer"] begin
		"--init_state"
    			help = "Initial seed for the random number generator (positive integer number)."
    			default = 45
			range_tester = check_is_uint64
    		"--init_seq"
    			help = "Identifier of the sequence produced by the "*
			    "random number generator (positive integer number)."
    			default = 54
			range_tester = check_is_uint64
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--num_of_rays"
			help = "number of `Ray`s generated for each integral evaluation"
			default = 10
			range_tester = check_is_uint64
		"--max_depth"
			help = "maximal number recursive integrations"
			default = 3
			range_tester = check_is_uint64
		"--russian_roulette_limit"
			help = "depth at whitch the Russian Roulette algorithm begins"
			default = 2
			range_tester = check_is_uint64
	end

	add_arg_group!(s["demo_animation"]["pointlight"], "pointlight renderer options");
	@add_arg_table! s["demo_animation"]["pointlight"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--ambient_color"
			help = "ambient color specified as '<R,G,B>' components. Example: --ambient_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--dark_parameter"
			help = "percentage of the retuned hit point color if it is not directly visible from a point-light source"
			default = 0.05
			range_tester = check_is_positive
	end

	#### RENDER #########################################################################92

	@add_arg_table! s["render"] begin
		"scenefile"
			help = "path to the file describing the scene to be rendered."
			arg_type = String
			#range_tester = input -> (typeof(query(input))<:File{format"PFM"})
			required = true
	end

	add_arg_group!(s["render"], "render options");
	@add_arg_table! s["render"] begin
		"--camera_type"
			help = "option for the camera type:\n"*
	    				"ort -> Orthogonal camera, per -> Perspective camera"
          	arg_type = String
			default = "per"
			range_tester = input -> check_is_one_of(input, CAMERAS)
		"--camera_position"
          	help = "camera position in the scene as '[X,Y,Z]'"
          	arg_type = String
          	default = "[-1,0,0]"
          	range_tester = check_is_vector
		"--alpha"
			help = "angle of view around z-axis, in degrees"
			arg_type = Float64
			default = 0.
		"--width"
			help = "pixel number on the width of the resulting demo image."
			default = 640
			range_tester = check_is_even_uint64
		"--height"
			help = "pixel number on the height of the resulting demo image."
			default = 480
			range_tester = check_is_even_uint64
          "--normalization", "-a"
			help = "scaling factor for the normalization process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.18
		"--gamma", "-g"
			help = "gamma value for the tone mapping process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 1.27
           "--avg_lum"
			help = "average luminosity of the resulting pfm image; if =0 (dafault value)"*
                    "it will be calculated automatically with the built-in avg_lum function"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.0
     	"--samples_per_pixel"
			help = "Number of samples per pixel for the antialiasing algorithm\n"*
					"It must be an integer perfect square, i.e. 0,1,4,9,16,...\n"*
					"If =0 (default value), antialiasing does not occurs."
     		default = 0
			range_tester = check_is_square
          "--set_pfm_name"
			help = "name of the pfm file to be saved"
			nargs = '?'
			arg_type = String
			default = "scene.pfm"
			constant = "scene.pfm"
		"--set_png_name"
			help = "name of the png file to be saved"
			nargs = '?'
			arg_type = String
			default = "scene.png"
			constant = "scene.png"
		"--ONLY_FOR_TESTS"
			help = "only for testing flag; do not use it"
			action = :store_true
	end

	add_arg_group!(s["render"],
		"render declare options for a scene; this options"*
		"allows to modify values of the scene without changing"*
		"the scenefile directly."
	);
	@add_arg_table! s["render"] begin
		"--declare_float"
			help = "Declare a variable. The syntax is «--declare-float=VAR:VALUE». Example: --declare_float=clock:150"
          	arg_type = String
			default = ""
			range_tester = check_is_declare_float
		end

	add_arg_group!(s["render"], "renderer to be used");
	@add_arg_table! s["render"] begin
		"onoff"
			action = :command
			help = "onoff renderer"
			#dest_name = "onoff"
		"flat"
			action = :command
			help = "flat-renderer"
			#dest_name = "flat"
		"pathtracer"
			action = :command
			help = "path tracing renderer"
			#dest_name = "pathtracer"
		"pointlight"
			action = :command
			help = "point-light tracing renderer"
			#dest_name = "pointlight"
	end

	add_arg_group!(s["render"]["onoff"], "onoff renderer options");
	@add_arg_table! s["render"]["onoff"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--color"
			help = "hit color specified as '<R,G,B>' components. Example: --ambient_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
	end

	add_arg_group!(s["render"]["flat"], "flat renderer options");
	@add_arg_table! s["render"]["flat"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
	end

	add_arg_group!(s["render"]["pathtracer"], "pathtracing renderer options");
	@add_arg_table! s["render"]["pathtracer"] begin
		"--init_state"
    			help = "Initial seed for the random number generator (positive integer number)."
    			default = 45
			range_tester = check_is_uint64
    		"--init_seq"
    			help = "Identifier of the sequence produced by the "*
			    "random number generator (positive integer number)."
    			default = 54
			range_tester = check_is_uint64
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--num_of_rays"
			help = "number of `Ray`s generated for each integral evaluation"
			default = 10
			range_tester = check_is_uint64
		"--max_depth"
			help = "maximal number recursive integrations"
			default = 3
			range_tester = check_is_uint64
		"--russian_roulette_limit"
			help = "depth at whitch the Russian Roulette algorithm begins"
			default = 2
			range_tester = check_is_uint64
	end

	add_arg_group!(s["render"]["pointlight"], "pointlight renderer options");
	@add_arg_table! s["render"]["pointlight"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--ambient_color"
			help = "ambient color specified as '<R,G,B>' components. Example: --ambient_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--dark_parameter"
			help = "percentage of the retuned hit point color if it is not directly visible from a point-light source"
			default = 0.05
			range_tester = check_is_positive
	end

	#### RENDER_ANIMATION #########################################################################92

	@add_arg_table! s["animation"] begin
		"--function"
			help = "name of the function that will be used to render the animation."
			arg_type = String
			range_tester = check_is_function
			required = true
		"--vec_variables"
			help = "vector of variable names that will change from frame to frame.\n"*
					"Must be declared as:  --vec_variables= \"[name1, name2, ...] \""
			arg_type = String
			range_tester = check_is_vec_variables
			required = true
		"--iterable"
			help = "iterable object from with the function will calcuate le variable values."
			arg_type = String
			range_tester = check_is_iterable
			required = true
		"scenefile"
			help = "path to the file describing the scene to be rendered."
			arg_type = String
			#range_tester = input -> (typeof(query(input))<:File{format"PFM"})
			required = true
	end

	add_arg_group!(s["animation"], "animation options");
	@add_arg_table! s["animation"] begin
		"--camera_type"
			help = "option for the camera type:\n"*
	    				"ort -> Orthogonal camera, per -> Perspective camera"
          	arg_type = String
			default = "per"
			range_tester = input -> check_is_one_of(input, CAMERAS)
		"--camera_position"
          	help = "camera position in the scene as '[X,Y,Z]'"
          	arg_type = String
          	default = "[-1,0,0]"
          	range_tester = check_is_vector
		"--alpha"
			help = "angle of view around z-axis, in degrees"
			arg_type = Float64
			default = 0.
		"--width"
			help = "pixel number on the width of the resulting demo image."
			default = 640
			range_tester = check_is_even_uint64
		"--height"
			help = "pixel number on the height of the resulting demo image."
			default = 480
			range_tester = check_is_even_uint64
          "--normalization", "-a"
			help = "scaling factor for the normalization process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.18
		"--gamma", "-g"
			help = "gamma value for the tone mapping process; must be positive"
			arg_type = Float64
			range_tester = check_is_positive
			default = 1.27
           "--avg_lum"
			help = "average luminosity of the resulting pfm image; if =0 (dafault value)"*
                    "it will be calculated automatically with the built-in avg_lum function"
			arg_type = Float64
			range_tester = check_is_positive
			default = 0.0
          "--set_anim_name"
			help = "name of the animation file to be saved"
			nargs = '?'
			arg_type = String
			default = "scene_animation.mp4"
			constant = "scene_animation.mp4"
     	"--samples_per_pixel"
			help = "Number of samples per pixel for the antialiasing algorithm\n"*
					"It must be an integer perfect square, i.e. 0,1,4,9,16,...\n"*
					"If =0 (default value), antialiasing does not occurs."
     		default = 0
			range_tester = check_is_square
		"--ONLY_FOR_TESTS"
			help = "only for testing flag; do not use it"
			action = :store_true
	end

	add_arg_group!(s["animation"],
		"animation declare options for a scene; this options"*
		"allows to modify values of the scene without changing"*
		"the scenefile directly."
	);
	@add_arg_table! s["animation"] begin
		"--declare_float"
			help = "Declare a variable. The syntax is «--declare-float=VAR:VALUE». Example: --declare_float=clock:150"
          	arg_type = String
			default = ""
			range_tester = check_is_declare_float
		end

	add_arg_group!(s["animation"], "renderer to be used");
	@add_arg_table! s["animation"] begin
		"onoff"
			action = :command
			help = "onoff renderer"
			#dest_name = "onoff"
		"flat"
			action = :command
			help = "flat-renderer"
			#dest_name = "flat"
		"pathtracer"
			action = :command
			help = "path tracing renderer"
			#dest_name = "pathtracer"
		"pointlight"
			action = :command
			help = "point-light tracing renderer"
			#dest_name = "pointlight"
	end

	add_arg_group!(s["animation"]["onoff"], "onoff renderer options");
	@add_arg_table! s["animation"]["onoff"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--color"
			help = "hit color specified as '<R,G,B>' components. Example: --ambient_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
	end

	add_arg_group!(s["animation"]["flat"], "flat renderer options");
	@add_arg_table! s["animation"]["flat"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
	end

	add_arg_group!(s["animation"]["pathtracer"], "pathtracing renderer options");
	@add_arg_table! s["animation"]["pathtracer"] begin
		"--init_state"
    			help = "Initial seed for the random number generator (positive integer number)."
    			default = 45
			range_tester = check_is_uint64
    		"--init_seq"
    			help = "Identifier of the sequence produced by the "*
			    "random number generator (positive integer number)."
    			default = 54
			range_tester = check_is_uint64
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--num_of_rays"
			help = "number of `Ray`s generated for each integral evaluation"
			default = 10
			range_tester = check_is_uint64
		"--max_depth"
			help = "maximal number recursive integrations"
			default = 3
			range_tester = check_is_uint64
		"--russian_roulette_limit"
			help = "depth at whitch the Russian Roulette algorithm begins"
			default = 2
			range_tester = check_is_uint64
	end

	add_arg_group!(s["animation"]["pointlight"], "pointlight renderer options");
	@add_arg_table! s["animation"]["pointlight"] begin
		"--background_color"
			help = "background color specified as '<R,G,B>' components. Example: --background_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--ambient_color"
			help = "ambient color specified as '<R,G,B>' components. Example: --ambient_color=<1,2,3>"
          	arg_type = String
          	default = "<0,0,0>"
          	range_tester = check_is_color
		"--dark_parameter"
			help = "percentage of the retuned hit point color if it is not directly visible from a point-light source"
			default = 0.05
			range_tester = check_is_positive
	end

	#### parse_args #####################################################################92

	parse_args(arguments, s)
end

function print_ArgParseSettings(vec::Dict{String,Any})
	parsed_command = vec["%COMMAND%"]
	parsed_args = vec[parsed_command]
     println("\nparsed_command : ", parsed_command)
	println("Parsed args:")
	for (key,val) in parsed_args
		println("  $key \t=>  $(repr(val))")
	end
	print("\n")
end


text:
## ExoplanetsSysSim/examples/run_abc.jl
## (c) 2015 Eric B. Ford

include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2017", "abc_setup_christiansen.jl"))
include(joinpath(Pkg.dir(),"ExoplanetsSysSim","examples","hsu_etal_2017", "process_abc_output.jl"))

using SysSimABC
using JLD

start_ind = 1
final_ind = 1

for n in start_ind:final_ind
  abc_plan = setup_abc(1)
  @time output = SysSimABC.run_abc(abc_plan)
  make_plots(output, EvalSysSimModel.get_ss_obs(), n)

  println("")
  println("Test ", n)
  println("")
  println("Mean = ")
  println(reshape(mean(output.theta,2), (1, size(output.theta,1))))
  println("")
  println("Standard deviation = ")
  println(reshape(std(output.theta,2), (1, size(output.theta,1))))
  println("")
  println(EvalSysSimModel.get_ss_obs())
  println("-----------------------------")
  save(string("test-",n,"-out.jld"), "output", output, "ss_true", EvalSysSimModel.get_ss_obs())
end

ext: jl


text:
@testset "Condition based iteration" begin
    @testset "Ambiguous nucleotide positions" begin
        function test_ns(A, seq)
            expected = Tuple{Int, eltype(A)}[]
            for i in 1:length(seq)
                if seq[i] == 'N'
                    push!(expected, (i, DNA_N))
                end
            end
            bioseq = LongSequence{A}(seq)
            @test collect(each(isambiguous, bioseq)) == expected
        end

        for len in [1, 10, 32, 1000, 10000, 100000]
            test_ns(DNAAlphabet{4}, random_dna(len))
            test_ns(RNAAlphabet{4}, random_rna(len))

            probs = [0.25, 0.25, 0.25, 0.25, 0.00]
            test_ns(DNAAlphabet{2}, random_dna(len, probs))
            test_ns(RNAAlphabet{2}, random_rna(len, probs))
        end

        dna_seq = dna"ANANANA"
        pos = [2, 4, 6]
        j = 1
        for (p, nuc) in each(isambiguous, dna_seq)
            @test p == pos[j]
            j += 1
        end

        dna_seq = dna"NATTCGRATY"
        pos = [1, 7, 10]
        j = 1
        for (p, nuc) in each(isambiguous, dna_seq)
            @test p == pos[j]
            j += 1
        end

        rna_seq = rna"ANANANA"
        pos = [2, 4, 6]
        j = 1
        for (p, nuc) in each(isambiguous, rna_seq)
            @test p == pos[j]
            j += 1
        end

        rna_seq = rna"NAUUCGRAUY"
        pos = [1, 7, 10]
        j = 1
        for (p, nuc) in each(isambiguous,rna_seq)
            @test p == pos[j]
            j += 1
        end
    end
end


text:
@testset "ImpulseResponse" begin
    df = exampledata(:gk)
    ns = (:logcpi, :logip, :ff, :ebp, :ff4_tc)
    r = lp(df, :ebp, wnames=ns, nlag=12, nhorz=48, vce=HRVCE())
    f = irf(r, :ebp, :ff4_tc, lag=1)

    @test sprint(show, f) == "ImpulseResponse"
    w = VERSION < v"1.4.0" ? " " : ""
    @test sprint(show, MIME("text/plain"), f, context=:displaysize=>(15,80)) == """
        ImpulseResponse with 48 horizons:
        ────────────────────────────────────────────────────────────────
             Estimate  Std. Error      z  Pr(>|z|)  Lower 90%  Upper 90%
        ────────────────────────────────────────────────────────────────
        0   0.394494     0.473835   0.83    0.4051  -0.384894   1.17388$w
        1   0.161872     0.509504   0.32    0.7507  -0.676188   0.999932
        2  -0.0203338    0.611236  -0.03    0.9735  -1.02573    0.98506$w
        3   0.152802     0.724653   0.21    0.8330  -1.03915    1.34475$w
        4   0.147288     0.737092   0.20    0.8416  -1.06512    1.3597$w$w
        5   1.0421       0.771451   1.35    0.1768  -0.226826   2.31102$w
        ────────────────────────────────────────────────────────────────
        42 rows of estimates omitted"""

    r1 = lp(df, :ebp, wnames=ns, nlag=12, minhorz=1, nhorz=5, vce=HRVCE())
    f1 = irf(r1, :ebp, :ff4_tc, lag=1)
    @test sprint(show, MIME("text/plain"), f1) == """
        ImpulseResponse with 5 horizons:
        ────────────────────────────────────────────────────────────────
             Estimate  Std. Error      z  Pr(>|z|)  Lower 90%  Upper 90%
        ────────────────────────────────────────────────────────────────
        1   0.161872     0.509504   0.32    0.7507  -0.676188   0.999932
        2  -0.0203338    0.611236  -0.03    0.9735  -1.02573    0.98506$w
        3   0.152802     0.724653   0.21    0.8330  -1.03915    1.34475$w
        4   0.147288     0.737092   0.20    0.8416  -1.06512    1.3597$w$w
        5   1.0421       0.771451   1.35    0.1768  -0.226826   2.31102$w
        ────────────────────────────────────────────────────────────────"""

    @test sprint(show, MIME("text/plain"), coeftable(f1, level=0.95)) == """
        ────────────────────────────────────────────────────────────────
             Estimate  Std. Error      z  Pr(>|z|)  Lower 95%  Upper 95%
        ────────────────────────────────────────────────────────────────
        1   0.161872     0.509504   0.32    0.7507  -0.836738    1.16048
        2  -0.0203338    0.611236  -0.03    0.9735  -1.21833     1.17767
        3   0.152802     0.724653   0.21    0.8330  -1.26749     1.5731$w
        4   0.147288     0.737092   0.20    0.8416  -1.29739     1.59196
        5   1.0421       0.771451   1.35    0.1768  -0.469918    2.55411
        ────────────────────────────────────────────────────────────────"""

    est = SmoothLP(:ff4_tc, search=grid(10), criterion=GCV())
    r2 = lp(est, df, :ebp, xnames=:ff4_tc, wnames=ns, nlag=12, minhorz=1, nhorz=5, vce=HARVCE(EWC()))
    f2 = irf(r2, :ebp, :ff4_tc)
    @test sprint(show, MIME("text/plain"), f2) == """
        ImpulseResponse with 5 horizons:
        ──────────────────────────────────────────────────────────────
           Estimate  Std. Error     z  Pr(>|z|)   Lower 90%  Upper 90%
        ──────────────────────────────────────────────────────────────
        1  0.657581    0.208515  3.15    0.0061   0.293537     1.02162
        2  0.870131    0.213599  4.07    0.0009   0.497212     1.24305
        3  0.948191    0.274859  3.45    0.0033   0.468319     1.42806
        4  0.892463    0.255986  3.49    0.0031   0.445541     1.33939
        5  0.703113    0.44229   1.59    0.1315  -0.0690743    1.4753$w
        ──────────────────────────────────────────────────────────────"""
end


text:
#
# Test script for least square solver
#
# author: Atsushi Sakai
#

using ArgMin

# ====Least square ====

function test_least_square()
	println("test_least_square")
	A = [ 0.97  1.86  0.41;
		 1.23  2.18  0.53;
		 0.80  1.24  0.62;
		 1.29  0.98  0.51;
		 1.10  1.23  0.69;
		 0.67  0.34  0.54;
		 0.87  0.26  0.62;
		 1.10  0.16  0.48;
		 1.92  0.22  0.71;
		 1.29  0.12  0.62];
	m, n = size(A)
	b = 1e3 * ones(m);
	xhat = solve_least_square(A, b)
	# println(xhat)

	@test xhat[1] ≈ 62.0766 atol=0.01
	@test xhat[2] ≈ 99.985 atol=0.01
	@test xhat[3] ≈ 1442.84 atol=0.01
end

function test_multi_objective_least_square()
	println("test_multi_objective_least_square")
	As = [[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0],
   		  [9.0, 2.0, 1.0, 4.0, 3.0, 1.6, 1.0, 2.0, -9.5, -19.0]];
	bs = [[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0],
   		  [5.0, 2.0, 8.0, -4.0, 1.0, 0.6, 17.0, 21.0, -9.0, -1.0]];
	lambdas = [0.1, 1.0]

	xhat = solve_multi_objective_least_square(As, bs, lambdas)

	@test xhat[1] ≈ 0.405977 atol=0.01
end

function test_solve_constrained_least_square()
	println("test_solve_constrained_least_square")
	A = randn(10,5)
	b = randn(10)
	C = randn(2,5)
	d = randn(2)
	xhat = solve_least_square(A,b,C,d);
end

function test_solve_nonlinear_least_square_with_newton_raphson()
	println("test_solve_nonlinear_least_square_with_newton_raphson")

	f(x) = (exp(x)-exp(-x)) / (exp(x)+exp(-x))
	Df(x) = 4.0 / (exp(x) + exp(-x))^2;

	xhat = solve_nonlinear_least_square_with_newton_raphson(f,Df,0.95)
end

function test_solve_nonlinear_least_square_with_gauss_newton()
	println("test_solve_nonlinear_least_square_with_gauss_newton")
	f(x) = (exp.(x) - exp.(-x)) ./ (exp.(x) + exp.(-x))
	Df(x) = 4 ./ (exp.(x) + exp.(-x)).^2
	xhat = solve_nonlinear_least_square_with_gauss_newton(f,Df,[1.15])
end

function test_solve_nonlinear_least_square_with_levenberg_marquardt()
	println("test_solve_nonlinear_least_square_with_levenberg_marquardt")
	f(x) = (exp.(x) - exp.(-x)) ./ (exp.(x) + exp.(-x))
	Df(x) = 4.0 ./ (exp.(x) + exp.(-x)).^2
	xhat = solve_nonlinear_least_square_with_levenberg_marquardt(f,Df,[0.95], 1.0)
end

function test_solve_constrained_nonlinear_least_square_with_augmented_lagragian()
	println("test_solve_constrained_nonlinear_least_square_with_augmented_lagragian")
	f(x) = [ x[1] + exp(-x[2]), x[1]^2 + 2*x[2] + 1 ]
	Df(x) = [ 1.0 -exp(-x[2]);
  			2*x[1]  2 ]
	g(x) = [ x[1] + x[1]^3 + x[2] + x[2]^2 ]
	Dg(x) = [ 1 + 3*x[1]^2  1 + 2*x[2] ]
	x, z, hist = solve_constrained_nonlinear_least_square_with_augmented_lagragian(f, Df, g, Dg, [0.5, -0.5], 1.0)
end

test_least_square()
test_multi_objective_least_square()
test_solve_constrained_least_square()
test_solve_nonlinear_least_square_with_newton_raphson()
test_solve_nonlinear_least_square_with_gauss_newton()
test_solve_nonlinear_least_square_with_levenberg_marquardt()
test_solve_constrained_nonlinear_least_square_with_augmented_lagragian()


text:
module Cronbach

using LinearAlgebra: diag

export cronbach

include("tauequivalent.jl")

end


text:
# Hyper parameters
ωs = [1.0, 0.0, 0.0]
ωp = 1e-7
# PO₄ mean and variance of observations fom WOA18
WOA = WorldOceanAtlasTools
μDIPobs3D, σ²DIPobs3D = WOA.fit_to_grid(grd, 2018, "phosphate", "annual", "1°", "an")
μDIPobs, σ²DIPobs = μDIPobs3D[iwet], σ²DIPobs3D[iwet]
μx = (μDIPobs, missing, missing)
σ²x = (σ²DIPobs, missing, missing)
δconvert(x) = @. exp(x)
cs = (δconvert, identity, identity)
# generate mismatch functions
f, ∇ₓf, ∇ₚf = generate_objective_and_derivatives(ωs, μx, σ²x, v, ωp)

#===========================================
Tests
===========================================#

@testset "Mismatch / cost functions" begin
    nt = length(T_all)
    nx = nt * nb
    np = length(p)
    x = p.xgeo * ones(nx)
    fval = f(x, p)
    ∇ₓfval = ∇ₓf(x, p)
    ∇ₚfval = ∇ₚf(x, p)
    @testset "Mismatch function" begin
        @test fval isa Float64
        @test fval > 0
    end
    @testset "Derivative w.r.t the state, x" begin
        @test size(∇ₓfval) == (1, nx)
    end
    @testset "Derivative w.r.t the parameters, p" begin
        @test size(∇ₚfval) == (1, np)
    end
end


text:
# %% Plot average winnings over time
using Distributions
using DataFrames
using Gadfly
using Colors
import Cairo, Fontconfig
using Printf
# %%
g = Geometric(0.5)

toss() = rand(g)+1

n_games = 20

n_players = 5
winnings = zeros(n_games,n_players)
for p in 1:n_players
    for n in 1:n_games
        payout = 0
        for i in 1:2^n
            payout += 2^toss()
        end
        winnings[n,p] = payout/(2^n)
    end
end

# %%

theme = Theme(major_label_font="PragmataPro Mono",
              minor_label_font="PragmataPro Mono",
              key_title_font="PragmataPro Mono Bold")
Gadfly.push_theme(theme)

x = 1:n_games
p = plot(x=x,
         Coord.cartesian(ymax=40),
         color=[RGBA(0,0,0,0.2)],
         Scale.x_continuous(labels = x -> @sprintf("2%d",x)),
         Guide.xlabel("ROUNDS PLAYED"),
         Guide.ylabel("MEAN PAYOUT"),
         Guide.title("Average Payout Per Round vs. Total Rounds Played"),
         )

for n in 1:n_players
    l = layer(y=winnings[:,n],
              Geom.line)
    push!(p,l)
end

draw(PDF("../../figures/average-over-time.pdf"), p)

legend = Guide.manual_color_key("", ["log2n+2"], ["deepskyblue"])
push!(p,legend)

l = layer(y=(1:n_games).+2,
          Geom.line,
          color=[theme.default_color],
          # color=["sine"],
          style(line_width=0.5mm),
          )

push!(p,l)
draw(PDF("../../figures/average-over-time-fit.pdf"), p)
# %%

max_stars_repo_name: akoen/petersburg

max_issues_repo_name: akoen/petersburg

max_forks_repo_name: akoen/petersburg


text:
"""
DecorrelationConfig contains config for the decorrelation steps
...
# Fields:
  - `steps::Int64`: number of decorrelation steps
  - `cov_estimate::Array{Float64,2}`: estimated (local) covariance,
  - `p_hr::Float64`: probability HR Steps
  - `p_mh::Float64`: probability MH Steps
...
"""
struct DecorrelationConfig
    steps::Int64
    cov_estimate::Array{Float64,2} # estimated (local) covariance,
    p_hr::Float64 # probability HR Steps
    p_mh::Float64 # probability MH Steps
end

"""
DecorrelationData
  collects data on performed decorrelations..
"""
struct DecorrelationData
    xv::Array{Float64,2} # all steps performed, chronologically (roughly)
    steps_success::Int64
    steps_fail::Int64
    steps::Vector{Int64} # +1 HR success, +2 MH success, -1 HR fail, -2 MH fail

    xx::Array{Float64,2} # on the way collected samples..
    xx_fdens::Array{Float64,1} # on the way collected samples..
    xx_fcost::Array{Float64,1} # on the way collected samples..

    sdata_hr::HRSamplingTrackingData
    sdata_mh::MHSamplingTrackingData
end

function DecorrelationData()
    return DecorrelationData( Array{Float64}(0,0) , Int64(0) , Int64(0) , Vector{Int64}() , Array{Float64}(0,0) , Array{Float64,1}() , Array{Float64,1}(), HRSamplingTrackingData() , MHSamplingTrackingData() )
end
function Base.hcat(a::DecorrelationData,b::DecorrelationData)
    return DecorrelationData( [a.xv b.xv] , a.steps_success+b.steps_success , a.steps_fail+b.steps_fail , [a.steps;b.steps] , [a.xx b.xx] , [a.xx_fdens ; b.xx_fdens] , [a.xx_fcost ; b.xx_fcost] , [ a.sdata_hr b.sdata_hr] , [ a.sdata_mh b.sdata_mh] )
end

#function DecorrelationData( dc1 , dc2 )
#    return DecorrelationData( [ dc1.xv dc2.xv ] , [dc1.steps_success ,dc2.steps_success ] , [dc1.steps_fail ,dc2.steps_fail ] , [dc1.xx ,dc2.xx ]  , [dc1.xxf ,dc2.xxf ]  )
#end

struct DecorrelationResult
    x::Vector{Float64}
    xfdens::Float64
    xfcost::Float64
    dd::DecorrelationData
end
function DecorrelationResult()
    return DecorrelationResult([NaN],NaN,NaN,DecorrelationData())
end

"""
CovEstimConfig configures the (local) covariance estimation
...
# Fields:
  - `n_considered::Int64`: the number of recently performed steps that we consider for estimating the covariance
  - `blowup_abs::Float64`: "Haario" regularization: add ( covestim_blowup_abs * diag(ones(d)) ) to the estimated cov. matrix
  - `blowup_rel_maxsv::Float64`: "Haario" rel regularization: add  ( covestim_blowup_abs * (diag(ones(d))) * norm(cov) ) , i.e. add relative to the frobenius norm ( max singular value)
  - `blowup_rel_minsv::Float64`: "Haario" rel regularization: add  ( covestim_blowup_abs * (diag(ones(d))) * min(svd(cov)[2]) ) , i.e. add relative to the smallest singular value
...

# Explanation about covestim_blowup parameters:
  blowup_abs is a standard procedure to regularize covariance estimation in adaptive mcmc. However, during
  ns, the estimated covariance will decrease, therefore, evetually, absolute blowup will prevent adaptiveness
  of the mcmc procedure.
  Therefore, we add the possibility to add blowup relative to singular values of the covariance matrix.
  The rel_maxsv blowup factor also limits the maximum non-isotropy of the estimated covariance in an intuitive way.
"""
struct CovEstimConfig
    n_considered::Int64 # the number of recently performed steps that we consider for estimating the covariance
    blowup_abs::Float64 # "Haario" regularization: add   ( covestim_blowup_abs * diag(ones(d)) ) to the estimated cov. matrix
    blowup_rel_maxsv::Float64 # "Haario" rel regularization: add  ( covestim_blowup_abs * (diag(ones(d))) * norm(cov) )        , i.e. add relative to the frobenius norm ( max singular value)
    blowup_rel_minsv::Float64 # "Haario" rel regularization: add  ( covestim_blowup_abs * (diag(ones(d))) * min(svd(cov)[2]) ) , i.e. add relative to the smallest singular value
end

# create a rather conservative estimation procedure by default
function CovEstimConfig(; n_considered=-1 , blowup_abs=NaN , blowup_rel_maxsv=NaN , blowup_rel_minsv=NaN)
    if(n_considered<0); n_considered = 200; end
    if(isnan(blowup_abs)); blowup_abs = 0.0001; end
    if(isnan(blowup_rel_maxsv)); blowup_rel_maxsv = 0.005; end # this also limits the ratio of small/large singular values to rouglhy 200
    if(isnan(blowup_rel_minsv)); blowup_rel_minsv = 0.1;   end # this also limits the ratio of small/large singular values
    return CovEstimConfig(n_considered,blowup_abs,blowup_rel_maxsv,blowup_rel_minsv)
end

"""
MCMCConfig configures the behavior of the mcmc decorrelation procedure
...
# Fields:
  - `p_hr::Float64`: probability of hr steps
  - `p_mh::Float64` : probability of mh steps
  - `cov_estim::CovEstimConfig` : covariance estimation configuration
  - `cov_initial::Array{Float64,2}` : initial covariance. Here, [NaN] is an alias for eye(dim)
...

"""
struct MCMCConfig
    p_hr::Float64
    p_mh::Float64
    cov_estim::CovEstimConfig
    cov_initial::Array{Float64,2} # if [NaN], then we take eye(dim)
end

function MCMCConfig()
    return MCMCConfig(0.5,0.5,CovEstimConfig(),reshape([NaN],1,1))
end

"""
NSConfig configures a run of the ns procedure
...
# Fields:
  - `steps_decorrelation::Int64`: the number of decorrelation steps
  - `shrinking_ratio::Float64` : the ratio by which we shrink in every iteration
  - `stop::StopCriterion` : the stop criterion
  - `mcmc::MCMCConfig` : configuration of the mcmc (decorrelation) procedure
...

"""
struct NSConfig
    steps_decorrelation::Int64
    shrinking_ratio::Float64
    stop::StopCriterion
    mcmc::MCMCConfig
    #decorrelation::DecorrelationConfig
    #save_live_samples::Bool
end

function NSConfig(; steps_decorrelation=-1 , shrinking_ratio=NaN , stop=StopCriterion() , mcmc=MCMCConfig() )#, save_live_samples=false)
    if(steps_decorrelation<0); steps_decorrelation=8; end
    if(isnan(shrinking_ratio)); shrinking_ratio=0.9; end
    return NSConfig(10,0.9,stop,mcmc)
end

struct NSIterData
    n_live_samples_a::Int64 # num samples before iteration
    n_live_samples_b::Int64 # num samples after removing sorted out samples
    n_live_samples_c::Int64 # num samples after replaing sorted out samples
    threshold::Float64      # cutoff threshold in this iteration (i.e. max fcost of samples "b" )
    X_sorted_out::Array{Float64} # sorted out samples
    Xf_sorted_out::Vector{Float64} # fvals of sorted out samples
    dd::Vector{DecorrelationData} # contains the data about the decorrelation steps
    ecov_raw::Array{Float64,2}
    ecov_regularized::Array{Float64,2}
end

struct NestedSamplingResult
    #X_sequence::Array{Float64,2}
    #X_sequence_f::Array{Float64,1}
    ns_iter_data::Array{NSIterData}
    #decorrelation_data::Array{DecorrelationData,1}
end

#function nested_sampling( X_init::Array{Float64,2} , f::Any , G::Array{Float64,2} , h::Array{Float64,1} , conf::NSConfig  ; X_f_init::Vector{Float64}=[] )

"""
  nested_sampling( X_init::Array{Float64,2} , f::Any , G::Array{Float64,2} , h::Array{Float64,1} , conf::NSConfig ; X_f_init::Vector{Float64}=[NaN] )

returns NestedSamplingResult

this function runs the nested sampling procedure
"""
function nested_sampling( X_init::Array{Float64,2} , f::Any , G::Array{Float64,2} , h::Array{Float64,1} , conf::NSConfig ; X_f_init::Vector{Float64}=[NaN] )

    output_level = 2

    X = copy(X_init)
    n_samples = size(X,2)
    d         = size(X,1)

    X_f::Array{Float64,1} = Array{Float64,1}()
    # init X_f if not given
    if(isnan(X_f_init[1]))
        X_f = Vector{Float64}(size(X_init,2))
        for zi in 1:size(X,2); X_f[zi] = f(X[:,zi]); end
    else
        X_f = copy(X_f_init)
    end

    # this data structure keeps track of the performed steps
    # create one for each step..
    #dd::Vector{DecorrelationData} = Vector{DecorrelationData}()

    itd::Vector{NSIterData}       = Vector{NSIterData}()

    # here we save the performed steps..
    Xv::Array{Float64,2}          = Array{Float64}(d,0)

    # this keeps track of all generated samples
    #X_sequence::Array{Float64,2}   = Array{Float64}(d,0)
    #X_sequence_f::Array{Float64,1} = Array{Float64}(0)

    # start the ns
    stop_ns = false
    t_start = time()
    iteration = 0
    while( !stop_ns )
        iteration += 1

        # profiling data
        iterdata_n_samples_a = size(X,2)

        # save live samples
        #X_saved_live_samples = Array{Float64,2}()
        #if(conf.save_live_samples)
        #    X_saved_live_samples   = X[:,:]
        #    X_saved_live_samples_f = X_f[:]
        #end

        # 1. sort out samples below threshold given by shrinking_ratio:
        println(X_f)
        f_cut = quantile( X_f , conf.shrinking_ratio)

        b_ok  = X_f .<= f_cut
        if(all(b_ok)); b_ok[ findmin(X_f)[2] ] = 0; end # at least sort out one.. :)

        #X_sequence   = [ X_sequence    X[ : , .~b_ok ] ]
        #X_sequence_f = [ X_sequence_f ; X_f[     .~b_ok ] ]
        iterdata_X_sorted_out   = X[ : , .~b_ok ]
        iterdata_Xf_sorted_out  = X_f[ .~b_ok ]

        #iterdata_fval_sorted_out = X_f[ .~b_ok ]

        X     = X[ : , b_ok ]
        X_f   = X_f[ b_ok ]

        iterdata_n_samples_b = size(X,2)

        # 2. select starting point samples
        n_replace = n_samples - size(X,2)
        idx_startingpoints = rand( 1:size(X,2) , n_replace )
        X0 = X[ : , idx_startingpoints ]

        # print iteration data..
        print( @sprintf("It %04d | %6.3f | %04d \n",iteration,f_cut,size(X,2)) )

        # 3. decorrelate
        # 3.1 estimate covariance..
        ecov_raw = eye(d)
        ecov     = NaN * ones(d,d)
        if( size(Xv,2)>4*d )
            if(size(Xv,2)>conf.mcmc.cov_estim.n_considered)
                ecov_raw = cov(Xv[ : , (end-conf.mcmc.cov_estim.n_considered+1):end ]' )
            else
                ecov_raw = cov(Xv')
            end
            sv_ecov = svd(ecov_raw)[2]
            ecov = ecov_raw + diagm(ones(d))*conf.mcmc.cov_estim.blowup_abs + diagm(ones(d)) * maximum(sv_ecov) * conf.mcmc.cov_estim.blowup_rel_maxsv + diagm(ones(d)) * minimum(sv_ecov) * conf.mcmc.cov_estim.blowup_rel_minsv

            if(output_level>=2)
                print("\necov_raw:\n$(ecov_raw)\n")
                print("ecov_regularized:\n$(ecov)\n")
                print("ecov_raw_sv: $(sv_ecov')\n")
                #print("ecov : blowup_abs: $(conf.mcmc.cov_estim.blowup_abs) , blowup_rel_max_sv: $(conf.mcmc.cov_estim.blowup_rel_maxsv*maximum(sv_ecov))) , blowup_rel_min_sv: $(conf.mcmc.cov_estim.blowup_rel_minsv*minimum(sv_ecov)) \n")
                print( @sprintf("ecov : blowup_abs: %6.3f  , blowup_rel_maxsv: %6.3f   , blowup_rel_minsv: %6.3f  \n", (conf.mcmc.cov_estim.blowup_abs) , (conf.mcmc.cov_estim.blowup_rel_maxsv*maximum(sv_ecov)) , (conf.mcmc.cov_estim.blowup_rel_minsv*minimum(sv_ecov)) ) )
            end
        else
            if( isnan(conf.mcmc.cov_initial[1]) )
                ecov = ecov_raw[:,:]
            else
                ecov = conf.mcmc.cov_initial[:,:]
            end
        end

        # create decorrelation data vector
        dd::Vector{DecorrelationData} = Vector{DecorrelationData}()

        threshold::Float64 = f_cut
        dc = DecorrelationConfig( conf.steps_decorrelation , ecov , conf.mcmc.p_hr , conf.mcmc.p_mh )
        sh_dcr::Array{DecorrelationResult,1} = Array{DecorrelationResult}(size(X0,2))
        # parallel / not parallel
        if(false)
            # 3.2 decorrelate samples..
            sh_dcr_af = Array{Future}(size(X0,2))
            for zi in 1:size(X0,2)
                sh_dcr_af[zi] = @spawn decorrelate( X0[:,zi] , f, threshold , G, h , dc )
            end
            for zi in 1:size(X0,2)
                sh_dcr[zi] = fetch(sh_dcr_af[zi])
            end
        else
            for zi in 1:size(X0,2)
                sh_dcr[zi] = run_mcmc( X0[:,zi] , NaN, Void(), f, threshold , G, h , dc )#decorrelate( X0[:,zi] , f, threshold , G, h , dc )
                #sh_dcr[zi] = decorrelate( X0[:,zi] , f, threshold , G, h , dc )
                #if(sh_dcr[zi].)
            end
        end

        X_new   = Array{Float64}(d,0)
        X_new_f = Array{Float64}(0)
        Xv_new  = Array{Float64}(d,0)
        for zi in 1:length(sh_dcr)
            push!(dd,sh_dcr[zi].dd)
            if(sh_dcr[zi].dd.steps_success>1)
                X_new   = [ X_new sh_dcr[zi].x ]
                Xv_new  = [ Xv_new  sh_dcr[zi].dd.xv ]
                X_new_f = [ X_new_f ; sh_dcr[zi].xfcost ]
            end
        end
        X   = [X     X_new   ]
        Xv  = [Xv    Xv_new  ]
        X_f = [X_f ; X_new_f ]

        iterdata_n_samples_c = size(X,2)

        push!(itd,NSIterData(iterdata_n_samples_a,iterdata_n_samples_b,iterdata_n_samples_c,f_cut,iterdata_X_sorted_out,iterdata_Xf_sorted_out,dd,ecov_raw,ecov))

        # check if we stop:
        if( time()-t_start>conf.stop.max_time)
            stop_ns = true

            # put remaining samples into the last iteration data
            push!(itd,NSIterData(iterdata_n_samples_c,0,0,-Inf,X,X_f,dd,ecov_raw,ecov))
        else
        end
    end

    #return NestedSamplingResult( [ X_sequence X ] , [ X_sequence_f;X_f ] , itd , dd )
    return NestedSamplingResult( itd  )
end

"""
  decorrelate( x_::Array{Float64,1} , f , threshold::Float64 , G::Array{Float64,2} , h::Array{Float64,1} , conf::DecorrelationConfig)

decorrelates the given samples inside {x|f(x)<threshold} by running mcmc steps

# Note
if all steps fail, then the returned value of x (xf) will be NaN

function decorrelate( x_::Array{Float64,1} , f , threshold::Float64 , G::Array{Float64,2} , h::Array{Float64,1} , conf::DecorrelationConfig)
    d = length(x_)

    dd_xx   = NaN * ones(d,conf.steps+1)
    dd_xx_f = NaN * ones(conf.steps+1)

    x  = copy(x_)
    xf = NaN

    hrconf = HRStepConfig()
    tdata_hr  = HRSamplingTrackingData()
    tdata_mh  = MHSamplingTrackingData()
    steps_data = Vector{Int64}(conf.steps)

    # config:
    mh_max_tries = 8

    cnt_successes = 0
    cnt_fails = 0
    rv_all = mvnrnd(x,conf.cov_estimate,conf.steps)
    step_types_all = sample_proportional([conf.p_hr;conf.p_mh],conf.steps)

    for zi in 1:conf.steps
        if(step_types_all[zi]==1)
            # HR STEP
            x_candidate_a = x - 4*rv_all[:,zi]  # times 4, because we want to be outside of viable space, if possible
            x_candidate_b = x + 4*rv_all[:,zi]
            #print("HRStep:\n")
            #print([x , x_candidate_a , x_candidate_b ])
            hr_result = hr_step( Void() , f , threshold , G , h , x , x_candidate_a, x_candidate_b , hrconf ,  tdata_hr  )
            if( hr_result[1]>0 )
                cnt_successes += 1
                steps_data[zi] = +1
                x  = hr_result[2]
                xf = hr_result[3]
                dd_xx[:,cnt_successes] = x
                dd_xx_f[cnt_successes] = xf
            else
                cnt_fails += 1
                steps_data[zi] = -1
            end
        end
        if(step_types_all[zi]==2)
            #MH STEP
            mh_result = mh_step( Void() , f , threshold , G , h , x , NaN , conf.cov_estimate , mh_max_tries , tdata_mh  )
            if(mh_result[1])
                cnt_successes += 1
                steps_data[zi] = +2
                x  = mh_result[2]
                xf = mh_result[4]
                dd_xx[:,cnt_successes] = x
                dd_xx_f[cnt_successes] = xf
            else
                cnt_fails += 1
                steps_data[zi] = -2
            end
        end
    end

    if(cnt_successes >= 1)
        # extract the successful steps data:
        dd_xx   = dd_xx[:,1:cnt_successes]
        dd_xx_f = dd_xx_f[1:cnt_successes]
        # extract xv data:
        xv = dd_xx - [ x_ dd_xx[:,1:(size(dd_xx,2)-1)] ]
        # create dc:
        dd = DecorrelationData(xv , cnt_successes , cnt_fails , steps_data , dd_xx,[NaN],dd_xx_f,tdata_hr,tdata_mh)
        return DecorrelationResult(x,NaN,xf,dd)
    else
        return DecorrelationResult(x,NaN,NaN,DecorrelationData( zeros(d,0) , cnt_successes , cnt_fails ,zeros(d,0),[NaN],tdata_hr,tdata_mh ))
    end
end
"""


text: using PkgBenchmark, GeometricFlux


text:
# A Frisch-Newton algorithm as described in http://projecteuclid.org/euclid.ss/1030037960
# qreg_ip_coef is a modified lp_fnm of Daniel Morillo & Roger Koenker
# Translated from Ox to Matlab by Paul Eilers 1999
# Modified by Roger Koenker 2000
# More changes by Paul Eilers 2004
# QuantileRegressions.jl developers 2015
#
# The functions in this file is released under the BSD-3 license, although the original
# versions were not. Permissing was given to Patrick Kofod Mogensen by Roger Koenker and
# Paul Eilers on April 24 2015 via E-mail. Original mail correspondences can be provided
# upon request.

function bound!(d, x, dx)
# Fill vector with allowed step lengths
# Replace with -x/dx for negative dx
    d .= 1e20
    @inbounds for i = 1:length(dx)
        if dx[i] < 0.0
            d[i] = -x[i] / dx[i]
        end
    end
    return d
end

function qreg_coef(Y, X::Matrix, p, method::IP)
# input: X is an n x k matrix of exogenous regressors,
#        Y is an n x 1 vector of outcome variables
#        p \in (0,1) is the quantile of interest
# Output: p^{th} regression quantiles.
# Construct the dual problem of quantile regression
c = -Y
T = eltype(X)

m = size(X, 1)
u = ones(m)
x = (1 - p) .* u
b = X'x

# Solve a linear program by the interior point method:
# min(c * u), s.t. A * x = b and 0 < x < u
# An initial feasible solution has to be provided as x

# Set some constants
beta = 0.99995
small = 1e-6
max_it = 50
n, m = size(X)
# Generate inital feasible point
s = u - x
y = -X\Y
dy = copy(y)
r = c - X*y
BLAS.axpy!(0.001, (r .== 0.0).*1.0, r)
z = r .* (r .> 0.0)
w = z .- r

gap = dot(c, x) - dot(y, b) + dot(w, u)

# set up caches
Xtmp = copy(X)
xinv, xi = copy(x), copy(x)
sinv = copy(s)
q = copy(z)
dx, ds, dz, dw = copy(w), copy(w), copy(w), copy(w)
fx, fs, fz, fw = copy(w), copy(w), copy(w), copy(w)
dxdz, dsdw = copy(w), copy(w)
tmp = copy(w)
# Start iterations
for it = 1:max_it
    #   Compute affine step
    @. q = 1 / (z / x + w / s)
    @. r = z - w
    Q = Diagonal(sqrt.(q)) # Very efficient to do since Q diagonal
    AQtF = qr(mul!(Xtmp, Q, X)) # PE 2004
    rhs = Q*r        # "
    ldiv!(dy, AQtF, rhs)
   # dy .= AQtF\rhs

    mul!(tmp, X, dy)
    @. dx = q*(tmp - r)
    @. ds = -dx
    @. dz = -z * (1 + dx / x)
    @. dw = -w * (1 + ds / s)

    # Compute maximum allowable step lengths
    bound!(fx, x, dx)
    bound!(fs, s, ds)
    bound!(fw, w, dw)
    bound!(fz, z, dz)
    @. tmp = min(fx, fs)
    fp = min(beta*minimum(tmp), 1)
    @. tmp = min(fw, fz)
    fd = min(beta*minimum(tmp), 1)
    # If full step is feasible, take it. Otherwise modify it
    if min(fp, fd) < 1.0

        # Update mu
        mu = dot(z, x) + dot(w, s)
        g = T(0)
        for i = 1:length(z)
          @inbounds g += (z[i] + fd[i]*dz[i])*(x[i] + fp[i]*dx[i]) + (w[i] + fd[i]*dw[i])*(s[i]+fp[i]*ds[i])
        end
        mu = mu * (g / mu)^3 / (2 * n)

        # Compute modified step
        @. dxdz = dx * dz
        @. dsdw = ds * dw
        @. xinv = 1 / x
        @. sinv = 1 / s
        @. xi = mu * (xinv - sinv)

        #rhs = rhs + Q * (dxdz - dsdw - xi)
        BLAS.axpy!(1.0, Q * (dxdz .- dsdw .- xi), rhs) # no gemv-wrapper gemv(Q, (dxdz - dsdw - xi), rhs,1,1,n)?
        ldiv!(dy, AQtF, rhs)

        mul!(tmp, X, dy)
        @. dx = q * (tmp + xi - r - dxdz + dsdw)
        @. ds = -dx
        for i = eachindex(dz, dw)
            @inbounds dz[i] = mu * xinv[i] - z[i] - xinv[i] * z[i] * dx[i] - dxdz[i]
            @inbounds dw[i] = mu * sinv[i] - w[i] - sinv[i] * w[i] * ds[i] - dsdw[i]
        end

        # Compute maximum allowable step lengths
        bound!(fx, x, dx)
        bound!(fs, s, ds)
        bound!(fw, w, dw)
        bound!(fz, z, dz)
        @. tmp = min(fx, fs)
        fp = min(beta*minimum(tmp), 1)
        @. tmp = min(fw, fz)
        fd = min(beta*minimum(tmp), 1)
    end

    # Take the steps
    BLAS.axpy!(fp, dx, x)
    BLAS.axpy!(fp, ds, s)
    BLAS.axpy!(fd, dy, y)
    BLAS.axpy!(fd, dw, w)
    BLAS.axpy!(fd, dz, z)

    gap = dot(c, x) - dot(y, b) + dot(w, u)

    if gap < small
        break
    end
end

return -y
end


text:
export Metric, AccuracyMetric, RankingMetric
export measure, count_true_positive

abstract type Metric end

abstract type AccuracyMetric <: Metric end

function measure(metric::AccuracyMetric, truth::AbstractVector, pred::AbstractVector)
    error("measure is not implemented for metric type $(typeof(metric))")
end

abstract type RankingMetric <: Metric end

function measure(metric::RankingMetric, truth::Array{T}, pred::Array{T}, k::Int) where T
    error("measure is not implemented for metric type $(typeof(metric))")
end

function count_true_positive(truth::Array{T}, pred::Array{T}) where T
    tp = 0
    for item in pred
        if findfirst(isequal(item), truth) != nothing
            tp += 1
        end
    end
    tp
end


text:
using Test
using Dates

using ScottishTaxBenefitModel

using .ModelHousehold:
    Household,
    Person,
    People_Dict,
    default_bu_allocation,
    get_benefit_units,
    get_head,
    get_spouse,
    is_single,
    num_people,
    num_children,
    pers_is_carer,
    pers_is_disabled,
    search

using .IncomeTaxCalculations:
    calc_income_tax!

using .NationalInsuranceCalculations:
    calculate_national_insurance!

using .Intermediate:
    MTIntermediate,
    apply_2_child_policy,
    make_intermediate

using .NonMeansTestedBenefits:
    calc_pre_tax_non_means_tested!,
    calc_post_tax_non_means_tested!

using .STBParameters:
    HoursLimits,
    HousingRestrictions,
    MinimumWage,
    UniversalCreditSys

using .UniversalCredit
    basic_conditions_satisfied,
    calc_elements!,
    calc_standard_allowance,
    calc_tariff_income!,
    calc_uc_child_costs!,
    calc_uc_income!,
    calc_universal_credit!,
    disqualified_on_capital,
    qualifiying_16_17_yo

using .Results:
    BenefitUnitResult,
    LMTResults,
    LMTCanApplyFor,
    init_household_result,
    init_benefit_unit_result,
    to_string

using .Utils:
    eq_nearest_p,
    to_md_table

using .GeneralTaxComponents:
    WEEKS_PER_MONTH,
    WEEKS_PER_YEAR

using .Definitions

using .STBIncomes
using .ExampleHelpers

## FIXME don't need both
sys = get_system( scotland=true )

@testset "UC Example Shakedown Tests" begin
    #
    # Just drive the example hhls through the UC routine
    # & see if anything crashes. Fuller tests to follow.
    #
    # Normally we'll do these tests monthly so they correspond better
    # to the CPAG examples
    # uc = get_default_uc( weekly = true )

    examples = get_all_examples()

    incomes = [110.0,145.0,325,755.0,1_000.0]

    for (hht,hh) in examples
        for income in incomes
            println( "on hhld '$hht' income=$income")

            bus = get_benefit_units( hh )
            intermed = make_intermediate(
                hh,
                sys.hours_limits,
                sys.age_limits,
                sys.child_limits )
            res = init_household_result( hh )
            hhead = get_head( hh )
            hhead.income[wages] = income
            calc_universal_credit!(
                res,
                hh,
                intermed,
                sys.uc,
                sys.age_limits,
                sys.hours_limits,
                sys.child_limits,
                sys.hr,
                sys.minwage
            )
            for buno in eachindex(res.bus)
                head = get_head( bus[buno] )
                println( res.bus[buno].uc )
                println( "UC Entitlement for $hht bu $buno earn $income = $(res.bus[buno].pers[head.pid].income[UNIVERSAL_CREDIT])" )
            end
        end
    end
end

@testset "Gail and Joe; CPAG 19/20 ch 3" begin
    # p41 example
    ucs = get_default_uc( weekly=false)
    g_and_j = get_example( cpl_w_2_children_hh )
    @test num_children( g_and_j ) == 2
    g_and_j.gross_rent = 0.0
    g_and_j.water_and_sewerage = 0.0
    g_and_j.mortgage_payment = 0.0
    g_and_j.other_housing_charges = 0.0
    g = get_head(g_and_j)
    j = get_spouse(g_and_j)

    set_childrens_ages!( g_and_j, 2, 3 ) # so none qualify for 1st child extra
    employ!(g)
    unemploy!(j)
    disable_seriously!(j)
    blind!(j)
    intermed = make_intermediate(
        g_and_j,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    @test intermed.buint[1].age_oldest_child == 3
    hres = init_household_result( g_and_j )
    calc_universal_credit!(
        hres,
        g_and_j,
        intermed,
        ucs,
        sys.age_limits,
        sys.hours_limits,
        sys.child_limits,
        sys.hr,
        sys.minwage
    )
    @test hres.bus[1].uc.standard_allowance == ucs.couple_oldest_25_plus
    @test hres.bus[1].uc.child_element ≈ 2*ucs.subsequent_child
    ## from j being disabled
    @test hres.bus[1].uc.limited_capacity_for_work_activity_element == ucs.limited_capcacity_for_work_activity
    # they have children, j is disabled, and no hcosts so..
    @test hres.bus[1].uc.work_allowance == ucs.work_allowance_no_housing
    println( hres.bus[1].uc )
    set_childrens_ages!( g_and_j, 10, 9 ) # so none qualify for 1st child extra
    intermed = make_intermediate(
        g_and_j,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    @test intermed.buint[1].age_oldest_child == 10
    calc_universal_credit!(
        hres,
        g_and_j,
        intermed,
        ucs,
        sys.age_limits,
        sys.hours_limits,
        sys.child_limits,
        sys.hr,
        sys.minwage
    )
    @test hres.bus[1].uc.child_element ≈ ucs.first_child+ucs.subsequent_child
    g_and_j.gross_rent = 100.0
    g_and_j.tenure = Private_Rented_Unfurnished
    intermed = make_intermediate(
        g_and_j,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    calc_universal_credit!(
        hres,
        g_and_j,
        intermed,
        ucs,
        sys.age_limits,
        sys.hours_limits,
        sys.child_limits,
        sys.hr,
        sys.minwage
    )
    println( "Tenure = $(g_and_j.tenure)" )
    @test hres.bus[1].uc.work_allowance == ucs.work_allowance_w_housing
    @test hres.bus[1].uc.housing_element ≈ 100.0
    # reset
    g_and_j.gross_rent = 0.0
    set_childrens_ages!( g_and_j, 3, 2 ) # so none qualify for 1st child extra
    intermed = make_intermediate(
        g_and_j,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    hres = init_household_result( g_and_j )
    hres.bus[1].pers[j.pid].income .= 0.0
    hres.bus[1].pers[g.pid].income .= 0.0
    empty!(g.assets)
    empty!(j.assets)
    g.assets[A_Savings_investments_etc] = 5_000.0
    hres.bus[1].pers[j.pid].income[CONTRIB_EMPLOYMENT_AND_SUPPORT_ALLOWANCE] = 483.81
    hres.bus[1].pers[g.pid].income[WAGES] = 1_200.0
    calc_universal_credit!(
        hres,
        g_and_j,
        intermed,
        ucs,
        sys.age_limits,
        sys.hours_limits,
        sys.child_limits,
        sys.hr,
        sys.minwage
    )
    println( hres.bus[1].uc )
    @test hres.bus[1].pers[j.pid].income[UNIVERSAL_CREDIT] ≈ 375.51
end

@testset "16-17yo adult tests; cpag ch3 sec2" begin
    ucs = get_default_uc( weekly=false)
    yph = get_example( mbu )
    head = get_head( yph )
    enable!( head )
    head.age = 17
    bus = get_benefit_units( yph )
    # doesn't qualify in standard case
    intermed = make_intermediate(
        yph,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    println( ModelHousehold.to_string( head ))
    println( Intermediate.to_string( intermed.buint[1] ))
        @test ! qualifiying_16_17_yo(
        bus[1],
        head,
        intermed.buint[1],
        ucs )
    # qualifies as carer
    carer!( head )
    intermed = make_intermediate(
        yph,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    @test qualifiying_16_17_yo(
        bus[1],
        head,
        intermed.buint[1],
        ucs )
    uncarer!( head )
    intermed = make_intermediate(
        yph,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    @test ! qualifiying_16_17_yo(
        bus[1],
        head,
        intermed.buint[1],
        ucs )

    add_child!( yph, 3, Female )
    bus = get_benefit_units( yph )
    @test num_people(bus[1]) == 2
    @test num_children(bus[1]) == 1
    intermed = make_intermediate(
        yph,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )
    @test qualifiying_16_17_yo(
        bus[1],
        head,
        intermed.buint[1],
        ucs )

end

@testset "capital" begin
    ucs = get_default_uc( weekly=false)
    cpl= get_example( cpl_w_2_children_hh )
    hres = init_household_result( cpl )
    head = get_head( cpl )
    spouse = get_spouse( cpl )
    empty!(head.assets)
    head.over_20_k_saving = false
    empty!(spouse.assets)
    spouse.over_20_k_saving = false
    bus = get_benefit_units( cpl )
    @test ! disqualified_on_capital(
        bus[1],
        ucs )
    spouse.over_20_k_saving = true
    @test disqualified_on_capital(
        bus[1],
        ucs )
    calc_tariff_income!(
        hres.bus[1],
        bus[1],
        ucs )
    @test hres.bus[1].uc.assets == 0
    @test hres.bus[1].uc.tariff_income == 0
    spouse.over_20_k_saving = false
    spouse.assets[A_Premium_bonds] = 8000
    @test ! disqualified_on_capital(
        bus[1],
        ucs )
        calc_tariff_income!(
            hres.bus[1],
            bus[1],
            ucs )
    @test hres.bus[1].uc.assets == 8_000
    @test hres.bus[1].uc.tariff_income ≈ ceil((8000-ucs.capital_min)/ucs.capital_tariff)
    head.assets[A_Premium_bonds] = 8_001
    @test disqualified_on_capital(
        bus[1],
        ucs )
    calc_tariff_income!(
        hres.bus[1],
        bus[1],
        ucs )
    @test hres.bus[1].uc.assets == 16_001
    @test hres.bus[1].uc.tariff_income ≈ ceil(ucs.capital_tariff\(16001-ucs.capital_min))

end

@testset "income calculations" begin
    ucs = get_default_uc( weekly=false)
    cpl= get_example( cpl_w_2_children_hh )
    head = get_head( cpl )
    head.age = 30
    spouse = get_spouse( cpl )
    spouse.age = 30
    empty!(head.assets)
    empty!(head.income)
    head.over_20_k_saving = false
    empty!(spouse.assets)
    empty!(spouse.income)
    spouse.over_20_k_saving = false
    bus = get_benefit_units( cpl )
    hres = init_household_result( cpl )
    intermed = make_intermediate(
        cpl,
        sys.hours_limits,
        sys.age_limits,
        sys.child_limits )

    calc_uc_income!(
        hres.bus[1],
        bus[1],
        intermed.buint[1],
        ucs,
        sys.minwage
    )
    @test hres.bus[1].uc.other_income == 0
    @test hres.bus[1].uc.earned_income == 0
    # CPAG 19/20 p120 example of low se income
    # I can't test this against the book exactly
    # because of weekly/monthly units but
    # this gets close.
    head.income[self_employment_income] = 300
    head.employment_status = Full_time_Self_Employed
    hres = init_household_result( cpl )
    @test hres.bus[1].pers[head.pid].income[SELF_EMPLOYMENT_INCOME] == 300
    hres.bus[1].pers[head.pid].income[INCOME_TAX] = 103.80
    minse = make_min_se(
        head.income[self_employment_income],
        head.age,
        ucs,
        sys.minwage
    )
    println( "minse=$minse")
    calc_uc_income!(
        hres.bus[1],
        bus[1],
        intermed.buint[1],
        ucs,
        sys.minwage
    )
    println( hres.bus[1].uc )
    @test hres.bus[1].uc.earned_income ≈ ucs.taper*(1249.9725-103.80-503.0)
end

run_full_tests = IS_LOCAL # && false

@testset "Run on actual Data" begin
    if run_full_tests
        #
        # Just a Shakedown run ..
        #
        nhhs,npeople = init_data()
        positive_uc = 0.0
        total_uc = 0.0
        for hno in 1:nhhs
            hh = get_household(hno)
            intermed = make_intermediate(
                hh,
                sys.hours_limits,
                sys.age_limits,
                sys.child_limits )

            hres = init_household_result( hh )
            println( "hhno $hno")
            # tax stuff, which we kinda sorta need
            bus = get_benefit_units( hh )
            calc_pre_tax_non_means_tested!(
                hres,
                hh,
                sys.nmt_bens,
                sys.hours_limits,
                sys.age_limits )

            for buno in eachindex(bus)
                # income tax, with some nonsense for
                # what remains of joint taxation..
                head = get_head( bus[buno] )
                spouse = get_spouse( bus[buno] )
                calc_income_tax!(
                    hres.bus[buno],
                    head,
                    spouse,
                    sys.it )
                for chno in bus[buno].children
                    child = bus[buno].people[chno]
                    calc_income_tax!(
                        hres.bus[buno].pers[child.pid],
                        child,
                        sys.it )
                end  # child loop
            end # bus loop
            calc_post_tax_non_means_tested!(
                hres,
                hh,
                sys.nmt_bens,
                sys.age_limits )

            calc_universal_credit!(
                hres,
                hh,
                intermed,
                sys.uc,
                sys.age_limits,
                sys.hours_limits,
                sys.child_limits,
                sys.hr,
                sys.minwage
            )
            for buno in eachindex(bus)
                # income tax, with some nonsense for
                # what remains of joint taxation..
                head = get_head( bus[buno] )
                uc = hres.bus[buno].pers[head.pid].income[UNIVERSAL_CREDIT]
                if uc > 0
                    positive_uc += hh.weight
                    total_uc += hh.weight*WEEKS_PER_YEAR*uc
                end
            end
        end # hhld loop
        total_uc /= 1_000_000
        println( "total UC entitlement $positive_uc total cost $total_uc (£m pa)")
        # mean of about 571pm ; actual 690.92 per StatXplore
    end # local
end # testset

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function setup_data()
    words = h5read(h5file, "s")
    worddict = Dict(words[i] => i for i=1:length(words))
    chardict = Dict{Char,Int}()

    traindoc = CoNLL.read(".data/wsj_00-18.conll")
    testdoc = CoNLL.read(".data/wsj_22-24.conll")
    info("# sentences of train doc: $(length(traindoc))")
    info("# sentences of test doc: $(length(testdoc))")

    traindata = setup_data(traindoc, worddict, chardict)
    testdata = setup_data(testdoc, worddict, chardict)
    info("# words: $(length(worddict))")
    info("# chars: $(length(chardict))")
    traindata, testdata, worddict, chardict
end

function setup_data(doc::Vector, worddict, chardict)
    data = []
    unkword = worddict["UNKNOWN"]
    root = worddict["PADDING"]
    for sent in doc
        w = Int[root]
        cs = Var[Var([0])]
        h = Int[0]
        for items in sent
            word = items[2]
            word0 = replace(word, r"[0-9]", '0')
            wordid = get(worddict, lowercase(word0), unkword)

            chars = Vector{Char}(word0)
            charids = map(c -> get!(chardict,c,length(chardict)+1), chars)

            head = parse(Int, items[7]) + 1
            push!(w, wordid)
            push!(cs, Var(charids))
            push!(h, head)
        end
        push!(data, (Var(w),cs,Var(h)))
    end
    data
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
export Interface
export handle

"""
An `Interface` belongs to a `FactorNode` and represents a half-edge.
An `Interface` has at most one partner interface, with wich it forms an edge.
"""
mutable struct Interface
    node::FactorNode
    edge::Union{AbstractEdge, Nothing}
    partner::Union{Interface, Nothing}
end

Interface(node::FactorNode) = Interface(node, nothing, nothing)

function show(io::IO, interface::Interface)
    iface_handle = handle(interface)
    (iface_handle == "") || (iface_handle = "($(iface_handle))")
    println(io, "Interface $(findfirst(isequal(interface), interface.node.interfaces)) $(iface_handle) of $(typeof(interface.node)) $(interface.node.id)")
end

"""Return interface handle name"""
function handle(interface::Interface)
    if isdefined(interface.node, :i)
        for h in keys(interface.node.i)
            if (typeof(h)==Symbol || typeof(h)==Int) && (interface.node.i[h] === interface)
                return string(h)
            end
        end
    end

    return ""
end

Base.isless(i1::Interface, i2::Interface) = isless(name(i1), name(i2))

name(iface::Interface) = string(iface.node.id)*handle(iface)
name(::Nothing) = ""

"""
Determines whether interface must be initialized with a breaker message;
i.e. for EP sites, loopy belief propagation, or situations where outbound
messages depend on inbound messages, such as a Nonlinear update without
a given inverse (RTS smoothing).
"""
function requiresBreaker(interface::Interface)
    partner_interface = ultimatePartner(interface)
    (partner_interface == nothing) && return false # Dangling edge

    return requiresBreaker(interface, partner_interface, partner_interface.node) # Dispatch to overloaded methods
end
requiresBreaker(interface::Interface, partner_interface::Interface, partner_node::FactorNode) = false # Default, function is overloaded for separate node types
requiresBreaker(::Nothing) = false # Failsafe

"""
Determine the type and dimensionality of the breaker interface message
"""
function breakerParameters(interface::Interface)
    partner_interface = ultimatePartner(interface)

    return breakerParameters(interface, partner_interface, partner_interface.node) # Dispatch to overloaded methods
end

"""
Constructs breaker types dictionary for breaker sites
"""
function breakerTypes(breaker_sites::Vector{Interface})
    breaker_types = Dict{Interface, Type}() # Initialize Interface to Message dictionary
    for site in breaker_sites
        (breaker_type, _) = breakerParameters(site)
        breaker_types[site] = breaker_type
    end

    return breaker_types
end


text:
function latexify(args...; kwargs...)

    ## Let potential recipes transform the arguments.
    args, kwargs = apply_recipe(args...; default_kwargs..., kwargs...)

    ## If the environment is unspecified, use auto inference.
    env = get(kwargs, :env, :auto)

    latex_function = infer_output(env, args...)

    result = latex_function(args...; kwargs...)

    COPY_TO_CLIPBOARD && clipboard(result)
    AUTO_DISPLAY && display(result)
    return result
end

apply_recipe(args...; kwargs...) = (args, kwargs)

# These functions should only be called from inside `latexify()`, so that
# `apply_recipe` gets a chance to change args
const OUTPUTFUNCTIONS = Dict(
                             :inline    => _latexinline,
                             :tabular   => _latextabular,
                             :table     => _latextabular,
                             :raw       => _latexraw,
                             :array     => _latexarray,
                             :align     => _latexalign,
                             :aligned   => (args...; kwargs...) -> _latexbracket(_latexalign(args...; kwargs..., aligned=true, starred=false); kwargs...),
                             :eq        => _latexequation,
                             :equation  => _latexequation,
                             :bracket   => _latexbracket,
                             :mdtable   => mdtable,
                             :mdtext    => mdtext,
                            )
function infer_output(env, args...)
    env === :auto && return get_latex_function(args...)
    # Must be like this, because items in OUTPUTFUNCTIONS must be defined
    env in [:arrows, :chem, :chemical, :arrow] && return _chemical_arrows
    return OUTPUTFUNCTIONS[env]
end

"""
    get_latex_function(args...)

Use overloading to determine which latex environment to output.

This determines the default behaviour of `latexify()` for different inputs.
"""
get_latex_function(args...) = _latexinline
get_latex_function(args::AbstractArray...) = _latexequation
get_latex_function(args::AbstractDict) = (args...; kwargs...) -> _latexequation(_latexarray(args...; kwargs...); kwargs...)
get_latex_function(args::Tuple...) = (args...; kwargs...) -> _latexequation(_latexarray(args...; kwargs...); kwargs...)
get_latex_function(arg::LaTeXString) = (arg; kwargs...) -> arg

function get_latex_function(x::AbstractArray{T}) where T <: AbstractArray
    try
        x = reduce(hcat, x)
        return (args...; kwargs...) -> _latexequation(_latexarray(args...; kwargs...); kwargs...)
    catch
        return _latexinline
    end
end

get_latex_function(lhs::AbstractVector, rhs::AbstractVector) = _latexalign


text:
struct NotMutable a end
x = NotMutable(2)
x.a = 3

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
struct Constraint{T}
    evaluate
    jacobian_state
    jacobian_action
    num_constraint::Int
    num_state::Int
    num_action::Int
    num_parameter::Int
    evaluate_cache::Vector{T}
    jacobian_state_cache::Matrix{T}
    jacobian_action_cache::Matrix{T}
    indices_inequality::Vector{Int}
end

Constraints{T} = Vector{Constraint{T}} where T

function Constraint(f::Function, num_state::Int, num_action::Int;
    indices_inequality::Vector{Int}=collect(1:0),
    num_parameter::Int=0)

    #TODO: option to load/save methods
    @variables x[1:num_state], u[1:num_action], w[1:num_parameter]

    evaluate = f(x, u, w)
    jacobian_state = Symbolics.jacobian(evaluate, x)
    jacobian_action = Symbolics.jacobian(evaluate, u)

    evaluate_func = eval(Symbolics.build_function(evaluate, x, u, w)[2])
    jacobian_state_func = eval(Symbolics.build_function(jacobian_state, x, u, w)[2])
    jacobian_action_func = eval(Symbolics.build_function(jacobian_action, x, u, w)[2])

    num_constraint = length(evaluate)

    return Constraint(
        evaluate_func,
        jacobian_state_func, jacobian_action_func,
        num_constraint, num_state, num_action, num_parameter,
        zeros(num_constraint), zeros(num_constraint, num_state), zeros(num_constraint, num_action),
        indices_inequality)
end

function Constraint()
    return Constraint(
        (c, x, u, w) -> nothing,
        (jx, x, u, w) -> nothing, (ju, x, u, w) -> nothing,
        0, 0, 0, 0,
        Float64[], Array{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0),
        collect(1:0))
end

function Constraint(f::Function, fx::Function, fu::Function, num_constraint::Int, num_state::Int, num_action::Int;
    indices_inequality::Vector{Int}=collect(1:0),
    num_parameter::Int=0)

    return Constraint(
        f,
        fx, fu,
        num_constraint, num_state, num_action, num_parameter,
        zeros(num_constraint), zeros(num_constraint, num_state), zeros(num_constraint, num_action),
        indices_inequality)
end

function constraint!(violations, constraints::Constraints{T}, states, actions, parameters) where T
    for (t, con) in enumerate(constraints)
        con.num_constraint == 0 && continue
        con.evaluate(con.evaluate_cache, states[t], actions[t], parameters[t])
        @views violations[t] .= con.evaluate_cache
        fill!(con.evaluate_cache, 0.0) # TODO: confirm this is necessary
    end
end

function jacobian!(jacobian_states, jacobian_actions, constraints::Constraints{T}, states, actions, parameters) where T
    H = length(constraints)
    for (t, con) in enumerate(constraints)
        con.num_constraint == 0 && continue
        con.jacobian_state(con.jacobian_state_cache, states[t], actions[t], parameters[t])
        @views jacobian_states[t] .= con.jacobian_state_cache
        fill!(con.jacobian_state_cache, 0.0) # TODO: confirm this is necessary
        t == H && continue
        con.jacobian_action(con.jacobian_action_cache, states[t], actions[t], parameters[t])
        @views jacobian_actions[t] .= con.jacobian_action_cache
        fill!(con.jacobian_action_cache, 0.0) # TODO: confirm this is necessary
    end
end


text:
struct IndicesZ
    q::Vector{Int}
    γ::Vector{Int}
    sγ::Vector{Int}
    ψ::Vector{Int}
    b::Vector{Int}
    sψ::Vector{Int}
    sb::Vector{Int}
end

function indices_z(model)
    nq = model.nq
    nc = model.nc
    q = collect(1:nq)
    γ = collect(nq .+ (1:nc))
    sγ = collect(nq + nc .+ (1:nc))
    ψ = collect(nq + nc + nc .+ (1:nc))
    b = collect(nq + nc + nc + nc .+ (1:nc))
    sψ = collect(nq + nc + nc + nc + nc .+ (1:nc))
    sb = collect(nq + nc + nc + nc + nc + nc .+ (1:nc))
    IndicesZ(q, γ, sγ, ψ, b, sψ, sb)
end

function initialize_z!(z, model, idx::IndicesZ, q)
    z[idx.q] .= q
    z[idx.γ] .= 1.0
    z[idx.sγ] .= 1.0
    z[idx.ψ] .= 1.0
    z[idx.b] .= 0.1
    z[idx.sψ] .= 1.0
    z[idx.sb] .= 0.1
end

struct Indicesθ
    q1::Vector{Int}
    q2::Vector{Int}
    u::Vector{Int}
    w::Vector{Int}
    f::Vector{Int}
    h::Vector{Int}
end

function indices_θ(model; nf=1)
    nq = model.nq
    nu = model.nu
    nw = model.nw

    q1 = collect(1:nq)
    q2 = collect(nq .+ (1:nq))
    u = collect(2nq .+ (1:nu))
    w = collect(2nq + nu .+ (1:nw))
    f = collect(2nq + nu + nw .+ (1:nf))
    h = collect(2nq + nu + nw + nf .+ (1:1))

    Indicesθ(q1, q2, u, w, f, h)
end

function initialize_θ!(θ, model, idx, q1, q2, u, w, f, h)
    θ[idx.q1] .= q1
    θ[idx.q2] .= q2
    θ[idx.u] .= u
    θ[idx.w] .= w
    θ[idx.f] .= f
    θ[idx.h] .= h
end

function indices_optimization(model)
    nq = model.nq
    nz = num_var(model)
    nc = model.nc
    IndicesOptimization(
        nz,
        nz,
        [collect(nq .+ (1:nc)), collect(nq + nc .+ (1:nc))],
        [collect(nq .+ (1:nc)), collect(nq + nc .+ (1:nc))],
        [[collect([nq + nc + nc + i, nq + nc + nc + nc + i]), collect([nq + nc + nc + nc + nc + i, nq + nc + nc + nc + nc + nc + i])] for i = 1:nc],
        [[collect([nq + nc + nc + i, nq + nc + nc + nc + i]), collect([nq + nc + nc + nc + nc + i, nq + nc + nc + nc + nc + nc + i])] for i = 1:nc],
        collect(1:(nq + nc + nc + nc)),
        collect(nq + nc + nc + nc .+ (1:nc)),
        collect(nq + nc + nc + nc + nc .+ (1:(2 * nc))),
        [collect(nq + nc + nc + nc + nc + (i - 1) * 2 .+ (1:2)) for i = 1:nc],
        collect(nq + nc + nc + nc .+ (1:(3 * nc))))
end


text:
"""
    min_even_theta_char(R)

Compute the even theta characteristic for which the theta constant at the input matrix in R has the smallest absolute value.
"""
function min_even_theta_char(R::RiemannMatrix)
    z = zeros(R.g);
    even_chars = even_theta_char(R.g);
    min_theta, char_index = findmin([abs(theta(z, R, char=c)) for c in even_chars])
    char = even_chars[char_index];
    return char
end

"""
    hessian(z, R, char=[])

Compute the hessian matrix of the theta function at input z and input matrix in R, with optional characteristics.
"""
function hessian(z::Array{<:Number}, R::RiemannMatrix, char::Array{}=[])
    M = zeros(Complex{Float64}, R.g, R.g);
    for i=1:R.g
        for j=i:R.g
            derivs = [zeros(R.g), zeros(R.g)];
            derivs[1][i] = 1;
            derivs[2][j] = 1;
            M[i,j] = theta(z, R, char=char, derivs=derivs);
            M[j,i] = M[i,j];
        end
    end
    return M
end

"""
    hessian(z, τ, char=[], siegel=true, ϵ=1.0e-12)

Compute the hessian matrix of the theta function at inputs z and τ, with optional characteristics, and optional inputs specifying whether to compute the Siegel transformation and the error in the value of the theta function.
"""
function hessian(z::Array{<:Number}, τ::Array{<:Number}, char::Array{}=[]; siegel::Bool=true, ϵ::Real=1.0e-12)
    R = RiemannMatrix(τ, siegel=siegel, ϵ=ϵ, nderivs=2);
    return hessian(z, R, char)
end

"""
    schottky_null(R, tol=1.0e-8)

Compute the even theta characteristic where the theta constant vanishes and the hessian matrix at the characteristic, as well as its nonzero eigenvalues, up to the input tolerance. Returns an array containing the characteristic, the hessian and the nonzero eigenvalues. If there is no such characteristic, return the absolute value of the theta constant at the smallest even characteristic.
"""
function schottky_null(R::RiemannMatrix, tol::Real=1.0e-8)
    char = min_even_theta_char(R);
    z = zeros(R.g);
    theta_value = abs(theta(z, R, char=char));
    if theta_value < tol # check that there is a vanishing theta null
        H = hessian(z, R, char);
        hessian_eigvals = filter(x -> abs(x) > tol, eigvals(H));
        return [char, H, hessian_eigvals]
    end
    return theta_value
end

"""
    schottky_null(τ, tol=1.0e-8)

Compute the even theta characteristic where the theta constant vanishes and the hessian matrix at the characteristic, as well as its nonzero eigenvalues, up to the input tolerance. Returns an array containing the characteristic, the hessian and the nonzero eigenvalues. If there is no such characteristic, return the absolute value of the theta constant at the smallest even characteristic.
"""
function schottky_null(τ::Array{<:Number}, tol::Real=1.0e-8)
    R = RiemannMatrix(τ, ϵ=tol, nderivs=2);
    return schottky_null(R, tol)
end


text:
function verify_storage_dir(folder::AbstractString = SERIALIZED_DIR)
    directory = abspath(normpath(folder))
    if !isdir(directory)
        mkpath(directory)
    end
end

function check_serialized_storage()
    verify_storage_dir(SERIALIZED_DIR)
    for path in SEARCH_DIRS
        verify_storage_dir(path)
    end
    return
end

function clear_serialized_system(name::String)
    seralized_file_extension =
        [".json", "_validation_descriptors.json", "_time_series_storage.h5"]
    file_names = [name * ext for ext in SERIALIZE_FILE_EXTENSIONS]
    for dir in SEARCH_DIRS
        for file in file_names
            if isfile(joinpath(dir, file))
                @debug "Deleting file" file
                rm(joinpath(dir, file), force = true)
            end
        end
    end
    return
end

function clear_all_serialized_system()
    for dir in SEARCH_DIRS
        @debug "Deleting dir" dir
        rm(dir, force = true, recursive = true)
    end
    check_serialized_storage()
    return
end

function get_serialization_dir(has_forecasts::Bool, has_reserves::Bool)
    if has_forecasts && has_reserves
        return SERIALIZE_FORECASTRESERVE_DIR
    elseif has_forecasts
        return SERIALIZE_FORECASTONLY_DIR
    elseif has_reserves
        return SERIALIZE_RESERVEONLY_DIR
    else
        return SERIALIZE_NOARGS_DIR
    end
end

get_serialized_filepath(name::String, has_forecasts::Bool, has_reserves::Bool) =
    joinpath(get_serialization_dir(has_forecasts, has_reserves), "$(name).json")

function is_serialized(name::String, has_forecasts::Bool, has_reserves::Bool)
    file_path = get_serialized_filepath(name, has_forecasts, has_reserves)
    if isfile(file_path)
        return true
    else
        return false
    end
end

function get_raw_data(; kwargs...)
    if haskey(kwargs, :raw_data)
        return kwargs[:raw_data]
    else
        throw(ArgumentError("Raw data directory not passed in build function."))
    end
end

function filter_kwargs(; kwargs...)
    system_kwargs = filter(x -> in(first(x), PSY.SYSTEM_KWARGS), kwargs)
    return system_kwargs
end

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
"""
    check_exited_exception(res::Exception) -> Bool

Recursively checks if an exception was caused by a worker exiting.
"""
check_exited_exception(res::CapturedException) =
    check_exited_exception(res.ex)
check_exited_exception(res::RemoteException) =
    check_exited_exception(res.captured)
check_exited_exception(res::ProcessExitedException) = true
check_exited_exception(res::Base.IOError) = true
check_exited_exception(res) = false

"""
    handle_fault(...)

An internal function to handle a worker dying or being killed by the OS.
Attempts to determine which `Thunk`s require rescheduling based on a
"deadlist", and then corrects the scheduler's internal `ComputeState` struct
to recover from the fault.

Note: The logic for this functionality is not currently perfectly robust to
all failure modes, and is only really intended as a last-ditch attempt to
repair and continue executing. While it should never cause incorrect execution
of DAGs, it *may* cause a `KeyError` or other failures in the scheduler due to
the complexity of getting the internal state back to a consistent and proper
state.
"""
function handle_fault(ctx, state, thunk, oldproc, chan, node_order)
    # Find thunks whose results were cached on the dead worker and place them
    # on what's called a "deadlist". This structure will direct the recovery
    # of the scheduler's state.
    deadlist = Thunk[thunk]
    # This thunk is guaranteed to not have valid cached data anymore
    thunk.cache = false
    thunk.cache_ref = nothing
    for t in keys(state.cache)
        v = state.cache[t]
        if v isa Chunk && v.handle isa DRef && v.handle.owner == oldproc.pid
            push!(deadlist, t)
            # Any inputs to dead cached thunks must be rescheduled
            function bfs!(deadlist, t)
                for input in t.inputs
                    istask(input) || continue
                    !(input in deadlist) && push!(deadlist, input)
                    bfs!(deadlist, input)
                end
            end
            bfs!(deadlist, t)
        end
    end
    # TODO: Find *all* thunks who were actively running on the dead worker

    # Empty cache of dead thunks
    for ct in keys(state.cache)
        if ct in deadlist
            delete!(state.cache, ct)
        end
    end

    function fix_waitdicts!(state, deadlist, t::Thunk; isleaf=false, offset=0)
        waiting, waiting_data = state.waiting, state.waiting_data
        off = repeat(" ", offset)
        offi = repeat(" ", offset+1)
        if !(t in keys(waiting))
            waiting[t] = Set{Thunk}()
        end
        if !isleaf
            # If we aren't a leaf thunk, then we may still need to recover
            # further into the DAG
            for input in t.inputs
                istask(input) || continue
                @assert haskey(waiting, t) "Error: $t not in state.waiting"
                push!(waiting[t], input)
                push!(waiting_data[input], t)
                isleaf = !(input in deadlist)
                fix_waitdicts!(state, deadlist, input; isleaf=isleaf, offset=offset+1)
            end
        end
        if isempty(waiting[t])
            delete!(waiting, t)
        end
    end

    # Add state.waiting deps back to state.waiting
    for ot in keys(state.waiting)
        fix_waitdicts!(state, deadlist, ot)
    end

    fix_waitdicts!(state, deadlist, thunk)

    # Remove thunks from state.ready that have inputs on the deadlist
    for idx in length(state.ready):-1:1
        rt = state.ready[idx]
        if any((input in deadlist) for input in rt.inputs)
            deleteat!(state.ready, idx)
        end
    end

    # Remove dead thunks from state.running, and add state.running
    # deps back to state.waiting
    wasrunning = copy(state.running)
    empty!(state.running)
    while !isempty(wasrunning)
        temp = pop!(wasrunning)
        if temp isa Thunk
            if !(temp in deadlist)
                push!(state.running, temp)
            end
            fix_waitdicts!(state, deadlist, temp)
        elseif temp isa Vector
            newtemp = []
            for t in temp
                fix_waitdicts!(state, deadlist, t)
                if !(t in deadlist)
                    push!(newtemp, t)
                end
            end
            isempty(newtemp) || push!(state.running, newtemp)
        else
            throw("Unexpected type in recovery: $temp")
        end
    end

    # Reschedule inputs from deadlist
    newproc = OSProc(rand(workers()))
    while length(deadlist) > 0
        dt = popfirst!(deadlist)
        if any((input in deadlist) for input in dt.inputs)
            # We need to schedule our input thunks first
            push!(deadlist, dt)
            continue
        end
        fire_task!(ctx, dt, newproc, state, chan, node_order)
        break
    end
end


text:
module GEXF

using GraphIO.EzXML
using Graphs
using Graphs: AbstractGraph, AbstractGraphFormat

import Graphs: savegraph

export GEXFFormat

# TODO: implement readgexf
struct GEXFFormat <: AbstractGraphFormat end
"""
    savegexf(f, g, gname)

Write a graph `g` with name `gname` to an IO stream `io` in the
[Gexf](http://gexf.net/format/) format. Return 1 (number of graphs written).
"""
function savegexf(io::IO, g::Graphs.AbstractGraph, gname::String)
    xdoc = XMLDocument()
    xroot = setroot!(xdoc, ElementNode("gexf"))
    xroot["xmlns"] = "http://www.gexf.net/1.2draft"
    xroot["version"] = "1.2"
    xroot["xmlns:xsi"] = "http://www.w3.org/2001/XMLSchema-instance"
    xroot["xsi:schemaLocation"] = "http://www.gexf.net/1.2draft/gexf.xsd"

    xmeta = addelement!(xroot, "meta")
    addelement!(xmeta, "description", gname)
    xg = addelement!(xroot, "graph")
    strdir = is_directed(g) ? "directed" : "undirected"
    xg["defaultedgetype"] = strdir

    xnodes = addelement!(xg, "nodes")
    for i in 1:nv(g)
        xv = addelement!(xnodes, "node")
        xv["id"] = "$(i-1)"
    end

    xedges = addelement!(xg, "edges")
    m = 0
    for e in Graphs.edges(g)
        xe = addelement!(xedges, "edge")
        xe["id"] = "$m"
        xe["source"] = "$(src(e)-1)"
        xe["target"] = "$(dst(e)-1)"
        m += 1
    end

    prettyprint(io, xdoc)
    return 1
end

savegraph(io::IO, g::AbstractGraph, gname::String, ::GEXFFormat) = savegexf(io, g, gname)

end #module


text:
primenumbers = [2]
num = 3
while true
  global num
  isprime = true
  for p in primenumbers
    if num % p == 0
      isprime = false
      break
    end
  end
  if isprime
    push!(primenumbers, num)
  end
  if length(primenumbers) > 10_000
    break
  end
  num += 2
end
println(maximum(primenumbers))

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl