text:
@testset "Crossover operators" begin

ss = symmetric_search_space(1, (0.0, 10.0))
fake_problem = FunctionBasedProblem(x -> 0.0, "test_problem", MinimizingFitnessScheme, ss)
DE = de_rand_1_bin(fake_problem, ParamsDict(
    :Population => reshape(collect(1.0:10.0), (1, 10)),
    :f => 0.4, :cr => 0.5, :NumParents => 3))

@testset "DiffEvoRandBin1" begin
    @testset "always copies from donor if length is 1" begin
        @test BlackBoxOptim.apply!(BlackBoxOptim.DiffEvoRandBin1(0.0, 0.0),
                                   [0.0], 4, DE.population, [1,2,3]) == [3.0]
    end

    @testset "always copies at least one element from donor" begin
        for i in 1:NumTestRepetitions
            len = rand(1:100)
            pop = rand(len,4)
            target = pop[:,1]
            saved_target = copy(target)
            res = BlackBoxOptim.apply!(BlackBoxOptim.DiffEvoRandBin1(0.0, 0.0),
                                       target, 1, pop, [2,3,4])
            @test (res === target)
            @test ndims( res ) == 1
            @test length( res ) == len
            @test sum(target .!= saved_target) == 1
        end
    end

    @testset "unlikely to copy everything if vectors are large" begin
        for i in 1:NumTestRepetitions
            len = 50 + rand(1:50)
            pop = rand(len,4)
            target = pop[:,1]
            saved_target = copy(target)
            res = BlackBoxOptim.apply!(BlackBoxOptim.DiffEvoRandBin1(0.1, 0.5),
                                       target, 1, pop, [2,3,4])
            @test any(target .!= saved_target)
            @test any(target .== saved_target)
        end
    end

    @testset "correctly modifies the parameters vector" begin
        @test BlackBoxOptim.apply!( BlackBoxOptim.DiffEvoRandBin1(1.0, 0.4),
                                                                [0.0], 4, DE.population, [1,2,3]) == [3.0 + (0.4 * (1.0 - 2.0))]
        @test BlackBoxOptim.apply!( BlackBoxOptim.DiffEvoRandBin1(1.0, 0.4),
                                                                [0.0], 4, DE.population, [2,3,1]) == [1.0 + (0.4 * (2.0 - 3.0))]
        @test BlackBoxOptim.apply!( BlackBoxOptim.DiffEvoRandBin1(1.0, 0.4),
                                                                [0.0], 4, DE.population, [3,2,1]) == [1.0 + (0.4 * (3.0 - 2.0))]
        @test BlackBoxOptim.apply!( BlackBoxOptim.DiffEvoRandBin1(1.0, 0.4),
                                                                [0.0], 4, DE.population, [1,3,5]) == [5.0 + (0.4 * (1.0 - 3.0))]
        @test BlackBoxOptim.apply!( BlackBoxOptim.DiffEvoRandBin1(1.0, 0.4),
                                                                [0.0], 3, DE.population, [4,9,8]) == [8.0 + (0.4 * (4.0 - 9.0))]

        pop2 = reshape(collect(1.0:8.0), 4, 2)'
        @test BlackBoxOptim.apply!( BlackBoxOptim.DiffEvoRandBin1(1.0, 0.6),
                                                                [0.0,0.0], 4, pop2, [1,2,3]) == [3.0 + (0.6 * (1.0 - 2.0)),
                                                                                7.0 + (0.6 * (5.0 - 6.0))]
        @test BlackBoxOptim.apply!( BlackBoxOptim.DiffEvoRandBin1(1.0, 0.6),
                                                                [0.0,0.0], 4, pop2, [1,2,4]) == [4.0 + (0.6 * (1.0 - 2.0)),
                                                                                8.0 + (0.6 * (5.0 - 6.0))]
    end
end

@testset "MutationWrapper" begin
    ss = symmetric_search_space(2, (0.0, 10.0))
    pop = reshape(collect(1.0:8.0), 2, 4)
    gibbs = UniformMutation(ss)
    gibbs_wrapper = BlackBoxOptim.MutationWrapper(gibbs)
    @test numchildren(gibbs_wrapper) == 1
    @test numparents(gibbs_wrapper) == 1
    mut_res = BlackBoxOptim.apply!(gibbs_wrapper, [0.0, 0.0], 1, pop, [2])
    @test sum(mut_res .== BlackBoxOptim.viewer(pop, 2)) == 0
end

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "distributions" begin
    include("test_distribution_functions.jl")
    include("test_bat_sampler.jl")
    include("test_HistogramAsUvDistribution.jl")
    include("test_HistogramAsMvDistribution.jl")
end

ext: jl


text:
using ContinuumArrays, Plots

###
# Dirichlet
####

L = LinearSpline(range(0,1; length=10_000))[:,2:end-1]
x = axes(L,1)
D = Derivative(x)
Δ = -((D*L)'D*L)
M = L'L
f = L \ exp.(x)
u = L * (Δ \ (M*f))
plot(u)

# we can also use cholesky
ũ = L * (cholesky(Symmetric(-Δ)) \ -(M*f))
plot(ũ)


text:
file_path = TEST_OUTPUTS

function test_plots(file_path::String; backend_pkg::String = "gr")
    if backend_pkg == "gr"
        Plots.gr()
    elseif backend_pkg == "plotlyjs"
        Plots.plotlyjs()
    else
        throw(error("$backend_pkg backend_pkg not supported"))
    end
    set_display = false
    cleanup = true
    @info("running tests with $backend_pkg with dispalay $set_display and cleanup $cleanup")

    (results_uc, results_ed) = run_test_sim(TEST_RESULT_DIR)
    problem_results = run_test_prob()
    gen_uc = PG.get_generation_data(results_uc)
    gen_ed = PG.get_generation_data(results_ed)
    gen_pb = PG.get_generation_data(problem_results)
    load_uc = PG.get_load_data(results_uc)
    load_ed = PG.get_load_data(results_ed)
    load_pb = PG.get_load_data(problem_results)
    svc_uc = PG.get_service_data(results_uc)
    svc_ed = PG.get_service_data(results_ed)
    svc_pb = PG.get_service_data(problem_results)

    @testset "test $backend_pkg plot production" begin
        out_path = joinpath(file_path, backend_pkg * "_plots")
        !isdir(out_path) && mkdir(out_path)
        plot_dataframe(
            gen_uc.data[:ActivePowerVariable__RenewableDispatch],
            gen_uc.time;
            set_display = set_display,
            title = "df_line",
            save = out_path,
        )
        plot_dataframe(
            gen_uc.data[:ActivePowerVariable__ThermalStandard],
            gen_uc.time,
            set_display = set_display,
            title = "df_stack",
            save = out_path,
            stack = true,
        )
        plot_dataframe(
            gen_uc.data[:ActivePowerVariable__ThermalStandard],
            gen_uc.time,
            set_display = set_display,
            title = "df_stair",
            save = out_path,
            stair = true,
        )
        plot_dataframe(
            gen_uc.data[:ActivePowerVariable__ThermalStandard],
            gen_uc.time,
            set_display = set_display,
            title = "df_bar",
            save = out_path,
            bar = true,
        )
        plot_dataframe(
            gen_uc.data[:ActivePowerVariable__ThermalStandard],
            gen_uc.time,
            set_display = set_display,
            title = "df_bar_stack",
            save = out_path,
            bar = true,
            stack = true,
        )
        plot_dataframe!(
            plot_dataframe(
                gen_uc.data[:ActivePowerVariable__ThermalStandard],
                gen_uc.time,
                set_display = set_display,
                stack = true,
            ),
            PG.no_datetime(load_uc.data[:Load]) .* -1,
            gen_uc.time,
            set_display = set_display,
            title = "df_gen_load",
            save = out_path,
        )

        list = readdir(out_path)
        expected_files = [
            "df_line.png",
            "df_stack.png",
            "df_stair.png",
            "df_bar.png",
            "df_bar_stack.png",
            "df_gen_load.png",
        ]
        # expected results not created
        @test isempty(setdiff(expected_files, list))
        # extra results created
        @test isempty(setdiff(list, expected_files))

        @info("removing test files")
        cleanup && rm(out_path, recursive = true)
    end

    @testset "test $backend_pkg pgdata plot production" begin
        out_path = joinpath(file_path, backend_pkg * "_pgdata_plots")
        !isdir(out_path) && mkdir(out_path)

        PG.plot_pgdata(
            gen_uc,
            set_display = set_display,
            title = "pg_data",
            save = out_path,
            bar = false,
            stack = false,
        )
        PG.plot_pgdata(
            gen_uc,
            set_display = set_display,
            title = "pg_data_stack",
            save = out_path,
            bar = false,
            stack = true,
        )
        PG.plot_pgdata(
            gen_uc,
            set_display = set_display,
            title = "pg_data_bar",
            save = out_path,
            bar = true,
            stack = false,
        )
        PG.plot_pgdata(
            gen_uc,
            set_display = set_display,
            title = "pg_data_bar_stack",
            save = out_path,
            bar = true,
            stack = true,
        )

        list = readdir(out_path)
        expected_files =
            ["pg_data.png", "pg_data_stack.png", "pg_data_bar.png", "pg_data_bar_stack.png"]
        # expected results not created
        @test isempty(setdiff(expected_files, list))
        # extra results created
        @test isempty(setdiff(list, expected_files))

        @info("removing test files")
        cleanup && rm(out_path, recursive = true)
    end

    @testset "test $backend_pkg demand plot production" begin
        out_path = joinpath(file_path, backend_pkg * "_demand_plots")
        !isdir(out_path) && mkdir(out_path)
        PG.plot_demand(
            results_uc,
            set_display = set_display,
            title = "demand",
            save = out_path,
            bar = false,
            stack = false,
            nofill = false,
        )
        PG.plot_demand(
            results_uc,
            set_display = set_display,
            title = "demand_stack",
            save = out_path,
            bar = false,
            stack = true,
            nofill = false,
        )
        PG.plot_demand(
            results_uc,
            set_display = set_display,
            title = "demand_bar",
            save = out_path,
            bar = true,
            stack = false,
            nofill = false,
        )
        PG.plot_demand(
            results_uc,
            set_display = set_display,
            title = "demand_bar_stack",
            save = out_path,
            bar = true,
            stack = true,
            nofill = false,
        )
        PG.plot_demand(
            results_uc,
            set_display = set_display,
            title = "demand_nofill",
            save = out_path,
            bar = false,
            stack = false,
            nofill = true,
        )
        PG.plot_demand(
            results_uc,
            set_display = set_display,
            title = "demand_nofill_stack",
            save = out_path,
            bar = false,
            stack = true,
            nofill = true,
        )
        PG.plot_demand(
            results_uc,
            set_display = set_display,
            title = "demand_nofill_bar",
            save = out_path,
            bar = true,
            stack = false,
            nofill = true,
        )
        PG.plot_demand(
            results_uc,
            set_display = set_display,
            title = "demand_nofill_bar_stack",
            save = out_path,
            bar = true,
            stack = true,
            nofill = true,
        )

        p = PG.plot_demand(
            results_uc.system,
            set_display = set_display,
            title = "sysdemand",
            save = out_path,
            aggregation = System,
        )
        plot_length = backend_pkg == "gr" ? length(p.series_list) : length(p.data)
        @test plot_length == 1

        p = PG.plot_demand(
            results_uc.system,
            set_display = set_display,
            title = "sysdemand_bus",
            save = out_path,
            aggregation = Bus,
        )
        plot_length = backend_pkg == "gr" ? length(p.series_list) : length(p.data)
        @test plot_length == 3

        list = readdir(out_path)
        expected_files = [
            "demand.png",
            "demand_stack.png",
            "demand_bar.png",
            "demand_bar_stack.png",
            "demand_nofill.png",
            "demand_nofill_stack.png",
            "demand_nofill_bar.png",
            "demand_nofill_bar_stack.png",
            "sysdemand.png",
            "sysdemand_bus.png",
        ]
        # expected results not created
        @test isempty(setdiff(expected_files, list))
        # extra results created
        @test isempty(setdiff(list, expected_files))

        @info("removing test files")
        cleanup && rm(out_path, recursive = true)
    end

    @testset "test $backend_pkg fuel plot production" begin
        out_path = joinpath(file_path, backend_pkg * "_fuel_plots")
        !isdir(out_path) && mkdir(out_path)

        PG.plot_fuel(
            results_uc,
            set_display = set_display,
            title = "fuel",
            save = out_path,
            bar = false,
            stack = false,
        )
        PG.plot_fuel(
            results_uc,
            set_display = set_display,
            title = "fuel_stack",
            save = out_path,
            bar = false,
            stack = true,
        )
        PG.plot_fuel(
            results_uc,
            set_display = set_display,
            title = "fuel_bar",
            save = out_path,
            bar = true,
            stack = false,
        )
        PG.plot_fuel(
            results_uc,
            set_display = set_display,
            title = "fuel_bar_stack",
            save = out_path,
            bar = true,
            stack = true,
        )

        list = readdir(out_path)
        expected_files =
            ["fuel.png", "fuel_stack.png", "fuel_bar.png", "fuel_bar_stack.png"]
        # expected results not created
        @test isempty(setdiff(expected_files, list))
        # extra results created
        @test isempty(setdiff(list, expected_files))

        @info("removing test files")
        cleanup && rm(out_path, recursive = true)
    end
end
try
    test_plots(file_path, backend_pkg = "gr")
    @info("done with GR, starting plotlyjs")
    test_plots(file_path, backend_pkg = "plotlyjs")
finally
    nothing
end

max_stars_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_issues_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_forks_repo_name: NREL-SIIP/InfrastructureSystemsPlots


text:
using ConfParser
using Test

conf = ConfParse(joinpath(@__DIR__, "confs", "config.http"))
parse_conf!(conf)

@test retrieve(conf, "email")    == "juliarocks@socks.com"
@test retrieve(conf, "password") == "qwerty"
@test commit!(conf, "email", "newemail@test.com") == true
@test haskey(conf, "email") == true
@test haskey(conf, "nothere") == false
save!(conf, joinpath(@__DIR__, "confs", "out.conf"))

conf = ConfParse(joinpath(@__DIR__, "confs", "out.conf"))
parse_conf!(conf)
@test retrieve(conf, "password") == "qwerty"
@test retrieve(conf, "email") == "newemail@test.com"


text:
module QGradMonomialBasesTests

using Test
using TensorPolynomialBases
using TensorValues

P = VectorValue{3}
V = VectorValue{3,Float64}
order = 1
basis = QGradMonomialBasis{P,V}(order)
G = gradient_type(basis)

x = VectorValue(2,3,5)
v = V[
  (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0),
  (3.0, 0.0, 0.0), (0.0, 5.0, 0.0), (0.0, 0.0, 2.0),
  (5.0, 0.0, 0.0), (0.0, 2.0, 0.0), (0.0, 0.0, 3.0),
  (15.0, 0.0, 0.0), (0.0, 10.0, 0.0), (0.0, 0.0, 6.0)]
g = G[
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0),
  (0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0),
  (0.0, 5.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 5.0, 0.0, 2.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0, 2.0, 0.0)]
test_polynomial_basis(basis,x,v,g)

P = VectorValue{2}
V = VectorValue{2,Float64}
order = 2
basis = QGradMonomialBasis{P,V}(order)
G = gradient_type(basis)

x = VectorValue(2.0,3.0)

v = V[
  (1.0, 0.0), (0.0, 1.0), (2.0, 0.0), (0.0, 3.0),
  (3.0, 0.0), (0.0, 2.0), (6.0, 0.0), (0.0, 6.0),
  (9.0, 0.0), (0.0, 4.0), (18.0, 0.0), (0.0, 12.0)]

g = G[
  (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (1.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 1.0), (0.0, 1.0, 0.0, 0.0), (0.0, 0.0, 1.0, 0.0),
  (3.0, 2.0, 0.0, 0.0), (0.0, 0.0, 3.0, 2.0), (0.0, 6.0, 0.0, 0.0),
  (0.0, 0.0, 4.0, 0.0), (9.0, 12.0, 0.0, 0.0), (0.0, 0.0, 12.0, 4.0)]
test_polynomial_basis(basis,x,v,g)

end # module


text:
const psxy  = plot
const psxy! = plot!
const psxyz  = plot3d
const psxyz! = plot3d!

# ---------------------------------------------------------------------------------------------------
function common_plot_xyz(cmd0::String, arg1, caller::String, first::Bool, is3D::Bool, kwargs...)
	arg2, arg3 = nothing, nothing
	N_args = (arg1 === nothing) ? 0 : 1
	is_ternary = (caller == "ternary") ? true : false
	if     (is3D)       gmt_proggy = (IamModern[1]) ? "plot3d "  : "psxyz "
	elseif (is_ternary) gmt_proggy = (IamModern[1]) ? "ternary " : "psternary "
	else		        gmt_proggy = (IamModern[1]) ? "plot "    : "psxy "
	end

	(occursin(" -", cmd0)) && return monolitic(gmt_proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode
	(!O) && (legend_type[1] = legend_bag())		# Make sure that we always start with an empty one

	cmd::String = "";	sub_module::String = ""	# Will change to "scatter", etc... if called by sub-modules
	g_bar_fill = Vector{String}()		# May hold a sequence of colors for gtroup Bar plots
	if (caller != "")
		if (occursin(" -", caller))		# some sub-modues use this piggy-backed call
			if ((ind = findfirst("|", caller)) !== nothing)	# A mixed case with "caler|partiall_command"
				sub_module = caller[1:ind[1]-1]
				cmd = caller[ind[1]+1:end]
				caller = sub_module		# Because of parse_BJR()
				(caller == "events") && (gmt_proggy = "events ")
			else
				cmd = caller
				caller = "others"		# It was piggy-backed
			end
		else
			sub_module = caller
			# Needs to be processed here to destinguish from the more general 'fill'
			(caller == "bar") && (g_bar_fill = helper_gbar_fill(d))
		end
	end

	if (occursin('3', caller) && !haskey(d, :p) && !haskey(d, :view) && !haskey(d, :perspective))
		d[:p] = "200/30"		# Need this before parse_BJR() so MAP_FRAME_AXES can be guessed.
	end

	if (is_ternary)
		cmd, opt_B::String = cmd * d[:B], d[:B]			# B option was parsed in plot/ternary
		delete!(d, :B)
		cmd, opt_R = parse_R(d, cmd, O)
	end
	if (is_ternary && !first) 	# Either a -J was set and we'll fish it here or no and we'll use the default.
		def_J = " -JX" * split(def_fig_size, '/')[1]
		cmd, opt_J::String = parse_J(d, cmd, def_J)
	else
		def_J = (is_ternary) ? " -JX" * split(def_fig_size, '/')[1] : ""		# Gives "-JX14c"
		if (is_ternary)  cmd, opt_J = parse_J(d, cmd, def_J)
		else             cmd, opt_B, opt_J, opt_R = parse_BJR(d, cmd, caller, O, def_J)
		end
	end

	cmd, opt_JZ = parse_JZ(d, cmd)
	cmd, = parse_common_opts(d, cmd, [:a :e :f :g :p :t :w :params], first)
	cmd, opt_l = parse_l(d, cmd)		# Parse this one (legend) aside so we can use it in classic mode
	cmd  = parse_these_opts(cmd, d, [[:D :shift :offset], [:I :intens], [:N :no_clip :noclip]])
	parse_ls_code!(d::Dict)				# Check for linestyle codes (must be before the GMTsyntax_opt() call)
	cmd  = GMTsyntax_opt(d, cmd)		# See if an hardcore GMT syntax string has been passed
	(is_ternary) && (cmd = add_opt(d, cmd, "M", [:M :dump]))
	opt_UVXY = parse_UVXY(d, "")		# Need it separate to not risk to double include it.
	cmd, opt_c = parse_c(d, cmd)		# Need opt_c because we may need to remove it from double calls

	# If a file name sent in, read it and compute a tight -R if this was not provided
	got_usr_R = (opt_R != "") ? true : false			# To know if the user set -R or we guessed it from data
	if (opt_R == "" && sub_module == "bar")  opt_R = "/-0.4/0.4/0"  end	# Make sure y_min = 0
	if (O && caller == "plotyy")
		cmd = replace(cmd, opt_R => "")					# Must remove old opt_R because a new one will be constructed
		ind = collect(findall("/", box_str[1])[2])		# 'box_str' was set in first call
		opt_R = '/' * box_str[1][4:ind[1]] * "?/?"		# Will become /x_min/x_max/?/?
	end
	cmd, arg1, opt_R, _, opt_i = read_data(d, cmd0, cmd, arg1, opt_R, is3D)
	(N_args == 0 && arg1 !== nothing) && (N_args = 1)
	(!O && caller == "plotyy") && (box_str[1] = opt_R)	# This needs modifications (in plotyy) by second call

	if (isGMTdataset(arg1) && getproj(arg1) == "" && opt_J == " -JX" * def_fig_size)
		cmd = replace(cmd, opt_J => " -JX" * split(def_fig_size, '/')[1] * "/0")	# If projected, it's a axis equal for sure
	end
	if (is3D && isempty(opt_JZ) && length(collect(eachmatch(r"/", opt_R))) == 5)
		cmd *= " -JZ6c"		# Default -JZ
	end

	cmd = add_opt(d, cmd, "A", [:A :steps :straight_lines], (x="x", y="y", meridian="m", parallel="p"))
	opt_F::String = add_opt(d, "", "", [:F :conn :connection],
	                (continuous=("c", nothing, 1), net=("n", nothing, 1), network=("n", nothing, 1), refpoint=("r", nothing, 1),  ignore_hdr="_a", single_group="_f", segments="_s", segments_reset="_r", anchor=("", arg2str)))
	(opt_F != "" && !occursin("/", opt_F)) && (opt_F = string(opt_F[1]))	# Allow con=:net or con=(1,2)
	(opt_F != "") && (cmd *= " -F" * opt_F)

	# Error Bars?
	got_Ebars = false
	val, symb = find_in_dict(d, [:E :error :error_bars], false)
	if (val !== nothing)
		cmd, arg1 = add_opt(add_opt, (d, cmd, "E", [symb]),
                            (x="|x",y="|y",xy="|xy",X="|X",Y="|Y", asym="_+a", colored="_+c", cline="_+cl", csymbol="_+cf", wiskers="|+n",cap="+w",pen=("+p",add_opt_pen)), false, isa(arg1, GMTdataset) ? arg1.data : (isa(arg1, Vector{<:GMTdataset}) ? arg1[1].data : arg1) )
		got_Ebars = true
		del_from_dict(d, [symb])
	end

	# Look for color request. Do it after error bars because they may add a column
	len = length(cmd);	n_prev = N_args;
	cmd, args, n, got_Zvars = add_opt(d, cmd, "Z", [:Z :level], :data, Any[arg1, arg2], (outline="_+l", fill="_+f"))
	if (n > 0)
		if (GMTver <= v"6.3")					# -Z is f again. Must save data into file to make it work.
			fname = joinpath(tempdir(), "GMTjl_temp_Z.txt")
			fid = open(fname, "w")
			for k = 1:length(args[n])  println(fid, args[n][k])  end;	close(fid)
			cmd *= fname
		else
			arg1, arg2 = args[:]
			N_args = n
		end
	end
	in_bag = (got_Zvars) ? true : false			# Other cases should add to this list
	if (N_args < 2)
		cmd, arg1, arg2, N_args = add_opt_cpt(d, cmd, CPTaliases, 'C', N_args, arg1, arg2, true, true, "", in_bag)
	else			# Here we know that both arg1 & arg2 are already occupied, so must use arg3 only
		cmd, arg3, = add_opt_cpt(d, cmd, CPTaliases, 'C', 0, arg3, nothing, true, true, "", in_bag)
		N_args = 3
	end

	mcc, bar_ok = false, (sub_module == "bar" && !check_bar_group(arg1))
	if ((!got_Zvars && !is_ternary) || bar_ok)	# If "bar" ONLY if not bar-group
		# See if we got a CPT. If yes there may be some work to do if no color column provided in input data.
		cmd, arg1, arg2, N_args, mcc = make_color_column(d, cmd, opt_i, len, N_args, n_prev, is3D, got_Ebars, bar_ok, g_bar_fill, arg1, arg2)
	end

	if (isempty(g_bar_fill))					# Otherwise bar fill colors are dealt somewhere else
		cmd = add_opt_fill(cmd, d, [:G :fill], 'G')
	end
	opt_Gsymb::String = add_opt_fill("", d, [:G :mc :markercolor :markerfacecolor :MarkerFaceColor], 'G')	# Filling of symbols

	# To track a still existing bug in sessions management at GMT lib level
	got_pattern = (occursin("-Gp", cmd) || occursin("-GP", cmd) || occursin("-Gp", opt_Gsymb) || occursin("-GP", opt_Gsymb)) ? true : false

	if (is_ternary)			# Means we are in the psternary mode
		cmd = add_opt(d, cmd, "L", [:L :vertex_labels])
	else
		opt_L::String = add_opt(d, "", "L", [:L :close :polygon],
		                (left="_+xl", right="_+xr", x0="+x", bot="_+yb", top="_+yt", y0="+y", sym="_+d", asym="_+D", envelope="_+b", pen=("+p",add_opt_pen)))
		(length(opt_L) > 3 && !occursin("-G", cmd) && !occursin("+p", cmd)) && (opt_L *= "+p0.5p")
		cmd *= opt_L
	end

	if ((val = find_in_dict(d, [:decorated])[1]) !== nothing)
		cmd = (isa(val, String)) ? cmd * " " * val : cmd * decorated(val)
	end

	opt_Wmarker::String = ""
	if ((val = find_in_dict(d, [:mec :markeredgecolor :MarkerEdgeColor])[1]) !== nothing)
		opt_Wmarker = "0.5p," * arg2str(val)		# 0.25p is too thin?
	end

	opt_W::String = add_opt_pen(d, [:W :pen], "W", true)     # TRUE to also seek (lw,lc,ls)
	((occursin("+c", opt_W)) && !occursin("-C", cmd)) &&
		@warn("Color lines (or fill) from a color scale was selected but no color scale provided. Expect ...")

	opt_S::String = add_opt(d, "", "S", [:S :symbol], (symb="1", size="", unit="1"))
	if (opt_S == "")			# OK, no symbol given via the -S option. So fish in aliases
		marca::String, arg1, more_cols = get_marker_name(d, arg1, [:marker, :Marker, :shape], is3D, true)
		if ((val = find_in_dict(d, [:ms :markersize :MarkerSize :size])[1]) !== nothing)
			(marca == "") && (marca = "c")			# If a marker name was not selected, defaults to circle
			if (isa(val, AbstractArray))
				(length(val) != size(arg1,1)) &&
					error("The size array must have the same number of elements rows in the data")
				arg1 = hcat(arg1, val[:])
			elseif (string(val) != "indata")
				marca *= arg2str(val);
			end
			opt_S = " -S" * marca
		elseif (marca != "")		# User only selected a marker name but no size.
			opt_S = " -S" * marca
			# If data comes from a file, then no automatic symbol size is added
			op = lowercase(marca[1])
			def_size = (op == 'p') ? "2p" : "7p"	# 'p' here stands for symbol points, not units
			(!more_cols && arg1 !== nothing && !isa(arg1, GMTcpt) && !occursin(op, "bekmrvw")) && (opt_S *= def_size)
		end
	else
		val, symb = find_in_dict(d, [:ms :markersize :MarkerSize :size])
		(val !== nothing) && @warn("option *$(symb)* is ignored when either *S* or *symbol* options are used")
		val, symb = find_in_dict(d, [:marker :Marker :shape])
		(val !== nothing) && @warn("option *$(symb)* is ignored when either *S* or *symbol* options are used")
	end

	opt_ML::String = ""
	if (opt_S != "")
		opt_ML, opt_Wmarker = parse_markerline(d, opt_ML, opt_Wmarker)
	end

	# See if any of the scatter, bar, lines, etc... was the caller and if yes, set sensible defaults.
	cmd  = check_caller(d, cmd, opt_S, opt_W, sub_module, g_bar_fill, O)
	(mcc && caller == "bar" && !got_usr_R && opt_R != " -R") && (cmd = recompute_R_4bars!(cmd, opt_R, arg1))	# Often needed
	_cmd = build_run_cmd(cmd, opt_B, opt_Gsymb, opt_ML, opt_S, opt_W, opt_Wmarker, opt_UVXY, opt_c)

	(got_Zvars && opt_S == "" && opt_W == "" && !occursin(" -G", _cmd[1])) && (_cmd[1] *= " -W0.5")

	# Let matrices with more data columns, and for which Color info was NOT set, plot multiple lines at once
	arg1 = helper_multi_cols(d, arg1, mcc, opt_R, opt_S, opt_W, caller, is3D, multi_col, _cmd,
	                         sub_module, g_bar_fill, got_Ebars, got_usr_R)

	# Try to limit the damage of this Fker bug in 6.2.0
	if ((mcc || got_Ebars) && (GMTver == v"6.2.0" && isGMTdataset(arg1) && occursin(" -i", cmd)) )
		if (isa(arg1, GMTdataset))	arg1 = arg1.data
		elseif (isa(arg1, Vector{<:GMTdataset}))
			(length(arg1) > 1) && @warn("Due to a bug in GMT6.2.0 I'm forced to use only the first segment")
			arg1 = arg1[1].data
		end
	end

	(!IamModern[1]) && put_in_legend_bag(d, _cmd, arg1, O, opt_l)

	_cmd = gmt_proggy .* _cmd				# In any case we need this
	_cmd = finish_PS_nested(d, _cmd)

	finish = (is_ternary && occursin(" -M",_cmd[1])) ? false : true		# But this case (-M) is bugged still in 6.2.0
	r = finish_PS_module(d, _cmd, "", K, O, finish, arg1, arg2, arg3)
	(got_pattern || occursin("-Sk", opt_S)) && gmt("destroy")  # Apparently patterns are screweing the session
	return r
end

# ---------------------------------------------------------------------------------------------------
function parse_markerline(d::Dict, opt_ML::String, opt_Wmarker::String)::Tuple{String, String}
	# Make this code into a function so that it can also be called from mk_styled_line!()
	if ((val = find_in_dict(d, [:ml :markerline :MarkerLine])[1]) !== nothing)
		if (isa(val, Tuple))           opt_ML = " -W" * parse_pen(val) # This can hold the pen, not extended atts
		elseif (isa(val, NamedTuple))  opt_ML = add_opt_pen(nt2dict(val), [:pen], "W")
		else                           opt_ML = " -W" * arg2str(val)
		end
		if (opt_Wmarker != "")
			@warn("markerline overrides markeredgecolor");		opt_Wmarker = ""
		end
	end
	return opt_ML, opt_Wmarker
end

# ---------------------------------------------------------------------------------------------------
function build_run_cmd(cmd, opt_B, opt_Gsymb, opt_ML, opt_S, opt_W, opt_Wmarker, opt_UVXY, opt_c)::Vector{String}
	# Build the executble command vector
	if (opt_W != "" && opt_S == "") 						# We have a line/polygon request
		_cmd = [cmd * opt_W * opt_UVXY]

	elseif (opt_W == "" && (opt_S != "" || opt_Gsymb != ""))	# We have a symbol request
		(opt_Wmarker != "" && opt_W == "") && (opt_Gsymb *= " -W" * opt_Wmarker)	# reuse var name
		(opt_ML != "") && (cmd *= opt_ML)					# If we have a symbol outline pen
		_cmd = [cmd * opt_S * opt_Gsymb * opt_UVXY]

	elseif (opt_W != "" && opt_S != "")						# We have both line/polygon and a symbol
		(occursin(opt_Gsymb, cmd)) && (opt_Gsymb = "")
		if (opt_S[4] == 'v' || opt_S[4] == 'V' || opt_S[4] == '=')
			_cmd = [cmd * opt_W * opt_S * opt_Gsymb * opt_UVXY]
		else
			(opt_Wmarker != "") && (opt_Wmarker = " -W" * opt_Wmarker)		# Set Symbol edge color
			cmd1 = cmd * opt_W * opt_UVXY
			(opt_B != " " && opt_B != "") && (cmd = replace(cmd, opt_B => ""))	# Some themes make opt_B = " "
			cmd2 = cmd * opt_S * opt_Gsymb * opt_Wmarker	# Don't repeat option -B
			(opt_c != "")  && (cmd2 = replace(cmd2, opt_c => ""))  				# Not in scond call (subplots)
			(opt_ML != "") && (cmd2 = cmd2 * opt_ML)				# If we have a symbol outline pen
			_cmd = [cmd1; cmd2]
		end

	else
		_cmd = [cmd * opt_UVXY]
	end
end

# ---------------------------------------------------------------------------------------------------
function helper_multi_cols(d::Dict, arg1, mcc, opt_R, opt_S, opt_W, caller, is3D, multi_col, _cmd, sub_module, g_bar_fill, got_Ebars, got_usr_R)
	# Let matrices with more data columns, and for which Color info was NOT set, plot multiple lines at once
	if (!mcc && opt_S == "" && (caller == "lines" || caller == "plot") && isa(arg1, Matrix{<:Real}) && size(arg1,2) > 2+is3D && size(arg1,1) > 1 && (multi_col[1] || haskey(d, :multicol)) )
		penC, penS = "", "";	cycle=:cycle;	multi_col[1] = false	# Reset because this is a use-only-once option
		(haskey(d, :multicol)) && delete!(d, :multicol)
		# But if we have a color in opt_W (idiotic) let it overrule the automatic color cycle in mat2ds()
		if (opt_W != "")  penT, penC, penS = break_pen(scan_opt(opt_W, "-W"))
		else              _cmd[1] *= " -W0.5"
		end
		if (penC != "")  cycle = [penC]  end
		arg1 = mat2ds(arg1, color=cycle, ls=penS, multi=true)	# Convert to multi-segment GMTdataset
		D::Vector{<:GMTdataset} = gmt("gmtinfo -C", arg1)		# But now also need to update the -R string
		_cmd[1] = replace(_cmd[1], opt_R => " -R" * arg2str(round_wesn(D[1].data)))
	elseif (!mcc && sub_module == "bar" && check_bar_group(arg1))	# !mcc because the bar-groups all have mcc = false
		_cmd[1], arg1 = bar_group(d, _cmd[1], opt_R, g_bar_fill, got_Ebars, got_usr_R, arg1)
	end
	return arg1
end

# ---------------------------------------------------------------------------------------------------
function helper_gbar_fill(d::Dict)::Vector{String}
	# This is a function that tryies to hammer the insistence that g_bar_fill is a Any
	# g_bar_fill may hold a sequence of colors for gtroup Bar plots
	gval = find_in_dict(d, [:fill :fillcolor], false)[1]	# Used for group colors
	if (isa(gval, Array{String}) && length(gval) > 1)
		g_bar_fill::Vector{String} = String[]
		append!(g_bar_fill, gval)
	elseif ((isa(gval, Array{Int}) || isa(gval, Tuple) && eltype(gval) == Int) && length(gval) > 1)
		g_bar_fill = Vector{String}(undef, length(gval))			# Patterns
		[g_bar_fill[k] = string('p',gval[k]) for k = 1:length(gval)]
	elseif (isa(gval, Tuple) && (eltype(gval) == String || eltype(gval) == Symbol) && length(gval) > 1)
		g_bar_fill = Vector{String}(undef, length(gval))			# Patterns
		[g_bar_fill[k] = string(gval[k]) for k = 1:length(gval)]
	else
		g_bar_fill = String[]		# To have somthing to return
	end
	return g_bar_fill
end

# ---------------------------------------------------------------------------------------------------
# Check if a group bar request or just bars. Returns TRUE in first case and FALSE in second
check_bar_group(arg1) = ( (isa(arg1, Matrix{<:Real}) || eltype(arg1) <: GMTdataset) &&
                          (isa(arg1, Vector{<:GMTdataset}) ? size(arg1[1],2) > 2 : size(arg1,2) > 2) )::Bool

# ---------------------------------------------------------------------------------------------------
function bar_group(d::Dict, cmd::String, opt_R::String, g_bar_fill::Array{String}, got_Ebars::Bool, got_usr_R::Bool, arg1)
	# Convert input array into a multi-segment Dataset where each segment is an element of a bar group
	# Example, plot two groups of 3 bars each: bar([0 1 2 3; 1 2 3 4], xlabel="BlaBla")

	if (got_Ebars)
		opt_E = scan_opt(cmd, "-E")
		((ind = findfirst("+", opt_E)) !== nothing) && (opt_E = opt_E[1:ind[1]-1])	# Strip eventual modifiers
		(((ind = findfirst("X", opt_E)) !== nothing) || ((ind = findfirst("Y", opt_E)) !== nothing)) && return cmd, arg1
		n_xy_bars = (findfirst("x", opt_E) !== nothing) + (findfirst("y", opt_E) !== nothing)
		n_cols = size(arg1,2)
		((n_cols - n_xy_bars) == 2) && return cmd, arg1			# Only one-bar groups
		(iseven(n_cols)) && error("Wrong number of columns in error bars array (or prog error)")
		n = Int((n_cols - 1) / 2)
		_arg = arg1[:, 1:(n+1)]				# No need to care with GMTdatasets because case was dealt in 'got_Ebars'
		bars_cols = arg1[:,(n + 2):end]		# We'll use this to appent to the multi-segments
	else
		_arg = isa(arg1, GMTdataset) ? deepcopy(arg1.data) : (isa(arg1, Vector{<:GMTdataset}) ? deepcopy(arg1[1].data) : deepcopy(arg1))
		bars_cols = missing
	end

	do_multi = true;	is_stack = false		# True for grouped; false for stacked groups
	is_hbar = occursin("-SB", cmd)				# An horizontal bar plot
	if ((val = find_in_dict(d, [:stack :stacked])[1]) !== nothing)
		# Take this (two groups of 3 bars) [0 1 2 3; 1 2 3 4]  and compute this (the same but stacked)
		# [0 1 0; 0 3 1; 0 6 3; 1 2 0; 1 5 2; 1 9 4]
		nl = size(_arg,2)-1				# N layers in stack
		tmp = zeros(size(_arg,1)*nl, 3)

		for m = 1:size(_arg, 1)			# Loop over number of groups
			tmp[(m-1)*nl+1,1] = _arg[m,1];		tmp[(m-1)*nl+1,2] = _arg[m,2];	# 3rd col is zero
			for n = 2:nl				# Loop over number of layers (n bars in a group)
				tmp[(m-1)*nl+n,1] = _arg[m,1]
				if (sign(tmp[(m-1)*nl+n-1,2]) == sign(_arg[m,n+1]))		# Because when we have neg & pos, case is diff
					tmp[(m-1)*nl+n,2] = tmp[(m-1)*nl+n-1,2] + _arg[m,n+1]
					tmp[(m-1)*nl+n,3] = tmp[(m-1)*nl+n-1,2]
				else
					tmp[(m-1)*nl+n,2] = _arg[m,n+1]
					tmp[(m-1)*nl+n,3] = 0
				end
			end
		end
		(is_hbar) && (tmp = [tmp[:,2] tmp[:,1] tmp[:,3]])	# Horizontal bars must swap 1-2 cols
		_arg = tmp
		do_multi = false;		is_stack = true
	end

	if (isempty(g_bar_fill) && findfirst("-G0/115/190", cmd) !== nothing)		# Remove the auto color
		cmd = replace(cmd, " -G0/115/190" => "")
	end

	# Convert to a multi-segment GMTdataset. There will be as many segments as elements in a group
	# and as many rows in a segment as the number of groups (number of bars if groups had only one bar)
	alpha = find_in_dict(d, [:alpha :fillalpha :transparency])[1]
	_argD = mat2ds(_arg; fill=g_bar_fill, multi=do_multi, fillalpha=alpha)
	(is_stack) && (_argD = ds2ds(_argD[1], fill=g_bar_fill, color_wrap=nl, fillalpha=alpha))
	if (is_hbar && !is_stack)					# Must swap first & second col
		for k = 1:length(_argD)
			_argD[k].data = [_argD[k].data[:,2] _argD[k].data[:,1]]
		end
	end
	(!isempty(g_bar_fill)) && delete!(d, :fill)

	if (bars_cols !== missing)		# Loop over number of bars in each group and append the error bar
		for k = 1:length(_argD)
			_argD[k].data = reshape(append!(_argD[k].data[:], bars_cols[:,k]), size(_argD[k].data,1), :)
		end
	end

	# Must fish-and-break-and-rebuild -S option
	opt_S = scan_opt(cmd, "-S")
	sub_b = ((ind = findfirst("+", opt_S)) !== nothing) ? opt_S[ind[1]:end] : ""	# The +Base modifier
	(sub_b != "") && (opt_S = opt_S[1:ind[1]-1])# Strip it because we need to (re)find Bar width
	bw = (isletter(opt_S[end])) ? parse(Float64, opt_S[3:end-1]) : parse(Float64, opt_S[2:end])	# Bar width
	n_in_group = length(_argD)					# Number of bars in the group
	gap = ((val = find_in_dict(d, [:gap])[1]) !== nothing) ? val/100 : 0	# Gap between bars in a group
	new_bw = (is_stack) ? bw : bw / n_in_group * (1 - gap)	# 'width' does not change in bar-stack
	new_opt_S = "-S" * opt_S[1] * "$(new_bw)u"
	cmd = (is_stack) ? replace(cmd, "-S"*opt_S*sub_b => new_opt_S*"+b") : replace(cmd, "-S"*opt_S => new_opt_S)

	if (!is_stack)								# 'Horizontal stack'
		g_shifts = linspace((-bw + new_bw)/2, (bw - new_bw)/2, n_in_group)
		col = (is_hbar) ? 2 : 1					# Horizontal and Vertical bars get shits in different columns
		for k = 1:n_in_group
			[_argD[k].data[n, col] += g_shifts[k] for n = 1:size(_argD[k].data,1)]
		end
	end

	if (!got_usr_R)								# Need to recompute -R
		info = gmt("gmtinfo -C", _argD)
		(info[1].data[3] > 0) && (info[1].data[3] = 0)		# If not negative then must be 0
		if (!is_hbar)
			dx = (info[1].data[2] - info[1].data[1]) * 0.005 + new_bw/2;
			dy = (info[1].data[4] - info[1].data[3]) * 0.005;
			info[1].data[1] -= dx;	info[1].data[2] += dx;	info[1].data[4] += dy;
			(info[1].data[3] != 0) && (info[1].data[3] -= dy);
		else
			dx = (info[1].data[2] - info[1].data[1]) * 0.005
			dy = (info[1].data[4] - info[1].data[3]) * 0.005 + new_bw/2;
			info[1].data[1] = 0.0;	info[1].data[2] += dx;	info[1].data[3] -= dy;	info[1].data[4] += dy;
			(info[1].data[1] != 0) && (info[1].data[1] -= dx);
		end
		info[1].data = round_wesn(info[1].data)		# Add a pad if not-tight
		new_opt_R = sprintf(" -R%.15g/%.15g/%.15g/%.15g", info[1].data[1], info[1].data[2], info[1].data[3], info[1].data[4])
		cmd = replace(cmd, opt_R => new_opt_R)
	end
	return cmd, _argD
end

# ---------------------------------------------------------------------------------------------------
function recompute_R_4bars!(cmd::String, opt_R::String, arg1)
	# Recompute the -R for bar plots (non-grouped), taking into account the width embeded in option S
	opt_S = scan_opt(cmd, "-S")
	sub_b = ((ind = findfirst("+", opt_S)) !== nothing) ? opt_S[ind[1]:end] : ""	# The +Base modifier
	(sub_b != "") && (opt_S = opt_S[1:ind[1]-1])# Strip it because we need to (re)find Bar width
	bw = (isletter(opt_S[end])) ? parse(Float64, opt_S[3:end-1]) : parse(Float64, opt_S[2:end])	# Bar width
	info = gmt("gmtinfo -C", arg1)
	dx::Float64 = (info[1].data[2] - info[1].data[1]) * 0.005 + bw/2;
	dy::Float64 = (info[1].data[4] - info[1].data[3]) * 0.005;
	info[1].data[1] -= dx;	info[1].data[2] += dx;	info[1].data[4] += dy;
	info[1].data = round_wesn(info[1].data)		# Add a pad if not-tight
	new_opt_R = sprintf(" -R%.15g/%.15g/%.15g/%.15g", info[1].data[1], info[1].data[2], 0, info[1].data[4])
	cmd = replace(cmd, opt_R => new_opt_R)
end

# ---------------------------------------------------------------------------------------------------
function make_color_column(d::Dict, cmd::String, opt_i::String, len::Int, N_args::Int, n_prev::Int, is3D::Bool, got_Ebars::Bool, bar_ok::Bool, bar_fill, arg1, arg2)
	# See if we got a CPT. If yes, there is quite some work to do if no color column provided in input data.
	# N_ARGS will be == n_prev+1 when a -Ccpt was used. Otherwise they are equal.

	(arg1 === nothing || isa(arg1, GMT.GMTcpt)) && return cmd, arg1, arg2, N_args, false  # Play safe

	mz, the_kw = find_in_dict(d, [:zcolor :markerz :mz])
	if ((!(N_args > n_prev || len < length(cmd)) && mz === nothing) && !bar_ok)	# No color request, so return right away
		return cmd, arg1, arg2, N_args, false
	end

	# Filled polygons with -Z don't need extra col
	((val = find_in_dict(d, [:G :fill], false)[1]) == "+z") && return cmd, arg1, arg2, N_args, false

	if     (isa(arg1, Vector{<:GMTdataset}))           n_rows, n_col = size(arg1[1])
	elseif (isa(arg1,GMTdataset) || isa(arg1, Array))  n_rows, n_col = size(arg1)
	end

	if (!isempty(bar_fill))
		if (isa(arg1,GMTdataset) || isa(arg1, Array))  arg1 = hcat(arg1, 1:n_rows)
		elseif (isa(arg1, Vector{<:GMTdataset}))       arg1[1].data = hcat(arg1[1].data, 1:n_rows)
		end
		arg2::GMTcpt = gmt(string("makecpt -T1/$(n_rows+1)/1 -C" * join(bar_fill, ",")))
		current_cpt[1] = arg2
		(!occursin(" -C", cmd)) && (cmd *= " -C")	# Need to inform that there is a cpt to use
		find_in_dict(d, [:G :fill])					# Must delete the :fill. Not used anymore
		return cmd, arg1, arg2, 2, true
	end

	warn1 = string("Probably color column in ", the_kw, " has incorrect dims. Ignoring it.")
	warn2 = "Plotting with color table requires adding one more column to the dataset but your -i
	option didn't do it, so you won't get what you expect. Try -i0-1,1 for 2D or -i0-2,2 for 3D plots"

	if (n_col <= 2+is3D)
		if (mz !== nothing)
			if (length(mz) != n_rows)  @warn(warn1); @goto noway  end
			if (isa(arg1,GMTdataset) || isa(arg1, Array))  arg1    = hcat(arg1, mz[:])
			elseif (isa(arg1, Vector{<:GMTdataset}))       arg1[1] = hcat(arg1[1], mz[:])
			end
		else
			if (opt_i != "")  @warn(warn2);		@goto noway		end
			cmd *= " -i0-$(1+is3D),$(1+is3D)"
			if ((val = find_in_dict(d, [:markersize :ms :size])[1]) !== nothing)
				cmd *= "-$(2+is3D)"		# Because we know that an extra col will be added later
			end
		end
	else
		if (mz !== nothing)				# Here we must insert the color col right after the coords
			if (length(mz) != n_rows)  @warn(warn1); @goto noway  end
			if (isa(arg1,GMTdataset) || isa(arg1, Array))  arg1    = hcat(arg1[:,1:2+is3D],    mz[:], arg1[:,3+is3D:end])
			elseif (isa(arg1, Vector{<:GMTdataset}))       arg1[1] = hcat(arg1[1][:,1:2+is3D], mz[:], arg1[1][:,3+is3D:end])
			end
		elseif (got_Ebars)				# The Error bars case is very multi. Don't try to guess then.
			if (opt_i != "")  @warn(warn2);	@goto noway  end
			cmd *= " -i0-$(1+is3D),$(1+is3D),$(2+is3D)-$(n_col-1)"
		end
	end

	if (N_args == n_prev)				# No cpt transmitted, so need to compute one
		if (mz !== nothing)                                    mi, ma = extrema(mz)
		else
			the_col = min(n_col,3)+is3D
			got_Ebars && (the_col -= 1)			# Bars => 2 cols
			if     (isa(arg1, Vector{<:GMTdataset}))           mi, ma = extrema(view(arg1[1], :, the_col))
			elseif (isa(arg1,GMTdataset) || isa(arg1, Array))  mi, ma = extrema(view(arg1,    :, the_col))
			end
		end
		just_C = cmd[len+2:end];	reset_i = ""
		if ((ind = findfirst(" -i", just_C)) !== nothing)
			reset_i = just_C[ind[1]:end]
			just_C  = just_C[1:ind[1]-1]
		end
		arg2 = gmt(string("makecpt -T", mi-0.001*abs(mi), '/', ma+0.001*abs(ma), " ", just_C) * (IamModern[1] ? " -H" : ""))
		current_cpt[1] = arg2
		if (occursin(" -C", cmd))  cmd = cmd[1:len+3]  end		# Strip the cpt name
		if (reset_i != "")  cmd *= reset_i  end		# Reset -i, in case it existed

		(!occursin(" -C", cmd)) && (cmd *= " -C")	# Need to inform that there is a cpt to use
		N_args = 2
	end

	@label noway

	return cmd, arg1, arg2, N_args, true
end

# ---------------------------------------------------------------------------------------------------
function get_marker_name(d::Dict, arg1, symbs::Vector{Symbol}, is3D::Bool, del::Bool=true)
	marca::String = "";		N = 0
	for symb in symbs
		if (haskey(d, symb))
			t = d[symb]
			if (isa(t, Tuple))				# e.g. marker=(:r, [2 3])
				msg = "";	cst = false
				o::String = string(t[1])
				if     (startswith(o, "E"))  opt = "E";  N = 3; cst = true
				elseif (startswith(o, "e"))  opt = "e";  N = 3
				elseif (o == "J" || startswith(o, "Rot"))  opt = "J";  N = 3; cst = true
				elseif (o == "j" || startswith(o, "rot"))  opt = "j";  N = 3
				elseif (o == "M" || startswith(o, "Mat"))  opt = "M";  N = 3
				elseif (o == "m" || startswith(o, "mat"))  opt = "m";  N = 3
				elseif (o == "R" || startswith(o, "Rec"))  opt = "R";  N = 3
				elseif (o == "r" || startswith(o, "rec"))  opt = "r";  N = 2
				elseif (o == "V" || startswith(o, "Vec"))  opt = "V";  N = 2
				elseif (o == "v" || startswith(o, "vec"))  opt = "v";  N = 2
				elseif (o == "w" || o == "pie" || o == "web" || o == "wedge")  opt = "w";  N = 2
				elseif (o == "W" || o == "Pie" || o == "Web" || o == "Wedge")  opt = "W";  N = 2
				end
				if (N > 0)  marca, arg1, msg = helper_markers(opt, t[2], arg1, N, cst)  end
				(msg != "") && error(msg)
				if (length(t) == 3 && isa(t[3], NamedTuple))
					if (marca == "w" || marca == "W")	# Ex (spiderweb): marker=(:pie, [...], (inner=1,))
						marca *= add_opt(t[3], (inner="/", arc="+a", radial="+r", size=("", arg2str, 1), pen=("+p", add_opt_pen)) )
					elseif (marca == "m" || marca == "M")
						marca *= vector_attrib(t[3])
					end
				end
			elseif (isa(t, NamedTuple))		# e.g. marker=(pie=true, inner=1, ...)
				key = keys(t)[1];	opt = ""
				if     (key == :w || key == :pie || key == :web || key == :wedge)  opt = "w"
				elseif (key == :W || key == :Pie || key == :Web || key == :Wedge)  opt = "W"
				elseif (key == :b || key == :bar)     opt = "b"
				elseif (key == :B || key == :HBar)    opt = "B"
				elseif (key == :l || key == :letter)  opt = "l"
				elseif (key == :K || key == :Custom)  opt = "K"
				elseif (key == :k || key == :custom)  opt = "k"
				elseif (key == :M || key == :Matang)  opt = "M"
				elseif (key == :m || key == :matang)  opt = "m"
				end
				if (opt == "w" || opt == "W")
					marca = opt * add_opt(t, (size=("", arg2str, 1), inner="/", arc="+a", radial="+r", pen=("+p", add_opt_pen)))
				elseif (opt == "b" || opt == "B")
					marca = opt * add_opt(t, (size=("", arg2str, 1), base="+b", Base="+B"))
				elseif (opt == "l")
					marca = opt * add_opt(t, (size=("", arg2str, 1), letter="+t", justify="+j", font=("+f", font)))
				elseif (opt == "m" || opt == "M")
					marca = opt * add_opt(t, (size=("", arg2str, 1), arrow=("", vector_attrib)))
				elseif (opt == "k" || opt == "K")
					marca = opt * add_opt(t, (custom="", size="/"))
				end
			else
				t1::String = string(t)
				(t1[1] != 'T') && (t1 = lowercase(t1))
				if     (t1 == "-" || t1 == "x-dash")    marca = "-"
				elseif (t1 == "+" || t1 == "plus")      marca = "+"
				elseif (t1 == "a" || t1 == "*" || t1 == "star")  marca = "a"
				elseif (t1 == "k" || t1 == "custom")    marca = "k"
				elseif (t1 == "x" || t1 == "cross")     marca = "x"
				elseif (is3D && (t1 == "u" || t1 == "cube"))  marca = "u"	# Must come before next line
				elseif (t1[1] == 'c')                   marca = "c"
				elseif (t1[1] == 'd')                   marca = "d"		# diamond
				elseif (t1 == "g" || t1 == "octagon")   marca = "g"
				elseif (t1[1] == 'h')                   marca = "h"		# hexagon
				elseif (t1 == "i" || t1 == "inverted_tri")  marca = "i"
				elseif (t1[1] == 'l')                   marca = "l"		# letter
				elseif (t1 == "n" || t1 == "pentagon")  marca = "n"
				elseif (t1 == "p" || t1 == "." || t1 == "point")  marca = "p"
				elseif (t1[1] == 's')                   marca = "s"		# square
				elseif (t1[1] == 't' || t1 == "^")      marca = "t"		# triangle
				elseif (t1[1] == 'T')                   marca = "T"		# Triangle
				elseif (t1[1] == 'y')                   marca = "y"		# y-dash
				elseif (t1[1] == 'f')                   marca = "f"		# for Faults in legend
				elseif (t1[1] == 'q')                   marca = "q"		# for Quoted in legend
				end
				t1 = string(t)		# Repeat conversion for the case it was lower-cased above
				# Still need to check the simpler forms of these
				if (marca == "")  marca = helper2_markers(t1, ["e", "ellipse"])   end
				if (marca == "")  marca = helper2_markers(t1, ["E", "Ellipse"])   end
				if (marca == "")  marca = helper2_markers(t1, ["j", "rotrect"])   end
				if (marca == "")  marca = helper2_markers(t1, ["J", "RotRect"])   end
				if (marca == "")  marca = helper2_markers(t1, ["m", "matangle"])  end
				if (marca == "")  marca = helper2_markers(t1, ["M", "Matangle"])  end
				if (marca == "")  marca = helper2_markers(t1, ["r", "rectangle"])   end
				if (marca == "")  marca = helper2_markers(t1, ["R", "RRectangle"])  end
				if (marca == "")  marca = helper2_markers(t1, ["v", "vector"])  end
				if (marca == "")  marca = helper2_markers(t1, ["V", "Vector"])  end
				if (marca == "")  marca = helper2_markers(t1, ["w", "pie", "web"])  end
				if (marca == "")  marca = helper2_markers(t1, ["W", "Pie", "Web"])  end
			end
			(del) && delete!(d, symb)
			break
		end
	end
	return marca, arg1, N > 0
end

function helper_markers(opt::String, ext, arg1, N::Int, cst::Bool)
	# Helper function to deal with the cases where one sends marker's extra columns via command
	# Example that will land and be processed here:  marker=(:Ellipse, [30 10 15])
	# N is the number of extra columns
	marca = "";	 msg = ""
	if (size(ext,2) == N && arg1 !== nothing)
		S = Symbol(opt)
		marca, arg1 = add_opt(add_opt, (Dict(S => (par=ext,)), opt, "", [S]), (par="|",), true, arg1)
	elseif (cst && length(ext) == 1)
		marca = opt * "-" * string(ext)
	else
		msg = string("Wrong number of extra columns for marker (", opt, "). Got ", size(ext,2), " but expected ", N)
	end
	return marca, arg1, msg
end

function helper2_markers(opt::String, alias::Vector{String})::String
	marca = ""
	if (opt == alias[1])			# User used only the one letter syntax
		marca = alias[1]
	else
		for k = 2:length(alias)		# Loop because of cases like ["w" "pie" "web"]
			o2 = alias[k][1:min(2,length(alias[k]))]	# check the first 2 chars and Ro, Rotrect or RotRec are all good
			if (startswith(opt, o2))  marca = alias[1]; break  end		# Good when, for example, marker=:Pie
		end
	end

	# If we still have found nothing, assume that OPT is a full GMT opt string (e.g. W/5+a30+r45+p2,red)
	(marca == "" && opt[1] == alias[1][1]) && (marca = opt)
	return marca
end

# ---------------------------------------------------------------------------------------------------
function check_caller(d::Dict, cmd::String, opt_S::String, opt_W::String, caller::String, g_bar_fill::Array{String}, O::Bool)::String
	# Set sensible defaults for the sub-modules "scatter" & "bar"
	if (caller == "scatter")
		if (opt_S == "")  cmd *= " -Sc5p"  end
	elseif (caller == "scatter3")
		if (opt_S == "")  cmd *= " -Su2p"  end
	elseif (caller == "lines")
		if (!occursin("+p", cmd) && opt_W == "") cmd *= " -W0.25p"  end # Do not leave without a pen specification
	elseif (caller == "bar")
		if (opt_S == "")
			bar_type = 0
			if (haskey(d, :bar))
				cmd, bar_opts = parse_bar_cmd(d, :bar, cmd, "Sb")
				bar_type = 1;	delete!(d, :bar)
			elseif (haskey(d, :hbar))
				cmd, bar_opts = parse_bar_cmd(d, :hbar, cmd, "SB")
				bar_type = 2;	delete!(d, :hbar)
			end
			if (bar_type == 0 || bar_opts == "")	# bar_opts == "" means only bar=true or hbar=true was used
				opt = (haskey(d, :width)) ? add_opt(d, "", "",  [:width]) : "0.8"	# The default
				_Stype = (bar_type == 2) ? " -SB" : " -Sb"
				cmd *= _Stype * opt * "u"

				optB = (haskey(d, :base)) ? add_opt(d, "", "",  [:base]) : "0"
				cmd *= "+b" * optB
			end
		end
		(isempty(g_bar_fill) && !occursin(" -G", cmd) && !occursin(" -C", cmd)) && (cmd *= " -G0/115/190")	# Default color
	elseif (caller == "bar3")
		if (haskey(d, :noshade) && occursin("-So", cmd))
			cmd = replace(cmd, "-So" => "-SO", count=1);
			delete!(d, :noshade)
		end
		if (!occursin(" -G", cmd) && !occursin(" -C", cmd))  cmd *= " -G0/115/190"	end
		if (!occursin(" -J", cmd))  cmd *= " -JX12c/0"  end
	end

	if (occursin('3', caller))
		if (!occursin(" -B", cmd) && !O)  cmd *= def_fig_axes3[1]  end	# For overlays default is no axes
	end

	return cmd
end

# ---------------------------------------------------------------------------------------------------
function parse_bar_cmd(d::Dict, key::Symbol, cmd::String, optS::String, no_u::Bool=false)::Tuple{String, String}
	# Deal with parsing the 'bar' & 'hbar' keywors of psxy. Also called by plot/bar3. For this
	# later module if input is not a string or NamedTuple the scatter options must be processed in bar3().
	# KEY is either :bar or :hbar
	# OPTS is either "Sb", "SB" or "So"
	# NO_U if true means to NOT automatic adding of flag 'u'
	opt::String = "";	got_str = false
	if (haskey(d, key))
		if (isa(d[key], String))
			opt, got_str = d[key], true
			cmd *= " -" * optS * opt;	delete!(d, key)
		elseif (isa(d[key], NamedTuple))
			opt = add_opt(d, "", optS, [key], (width="",unit="1",base="+b",height="+B",nbands="+z",Nbands="+Z"))
		elseif (isa(d[key], Bool) && d[key])
		else
			error("Argument of the *bar* keyword can be only a string or a NamedTuple.")
		end
	end

	if (opt != "" && !got_str)				# Still need to finish parsing this
		flag_u = no_u ? "" : 'u'
		if ((ind = findfirst("+", opt)) !== nothing)	# See if need to insert a 'u'
			if (!isletter(opt[ind[1]-1]))  opt = opt[1:ind[1]-1] * flag_u * opt[ind[1]:end]  end
		else
			pb = (optS != "So") ? "+b0" : ""		# The default for bar3 (So) is set in the bar3() fun
			if (!isletter(opt[end]))  opt *= flag_u	  end	# No base set so default to ...
			opt *= pb
		end
		cmd *= opt
	end
	return cmd, opt
end


text:
module BanditOpt

"""
BanditOpt is a Julia package for experimenting with various Multi-Armed Bandit algorithms.
"""
BanditOpt

"""
BanditOpt.Agents module provides implementations of various multi-arm bandit algorithms and an
unified API to access those algorithms. This includes base types for implementing various types
of agents and the common interface.
"""
module Agents
    using Distributions
    using LaTeXStrings
    using Printf
    include( "Agents/AgentBase.jl" )
    include( "Agents/StationaryAgentBase.jl" )

    export
        # Methods
        get_arm_index, update_reward!, reset!, info_str,
        show,
        # Agents
        AgentBase,
        # Stationary Agents
        StationaryAgentBase,
        epsGreedy, epsNGreedy
end

"""
    Arms

`Arms` module includes the different types of reward processes available for
Bandit algorithms. It define the following base types for creating new reward
processes.
- [Arms.ArmBase](@ref)
- [Arms.StationaryArmBase](@ref)
- [Arms.NonstationaryArmBase](@ref)
The types and methods in `Arms` are not exported and need to be called e.g as
`Arms.pull!()`
"""
module Arms
    using Distributions
    include( "Arms/ArmBase.jl" )
    include( "Arms/StationaryArms.jl" )
    include( "Arms/NonStationaryArms.jl" )

    # export
    #     # Methods
    #     pull!, tick!, reset!,
    #     # Arm Models
    #     ArmBase,
    #     # Stationary Arm Models
    #     StationaryArmBase,
    #     # Non-stationary Arm Models
    #     NonstationaryArmBase
end

# -----------------------------  Exports ----------------------------------
using .Agents
using .Arms

export
    # Agents Module
    Agents,
        # Methods
        get_arm_index, update_reward!, reset!, info_str, show,
        # Agents
        AgentBase,
        # Stationary Agents
        StationaryAgentBase,
        epsGreedy, epsNGreedy,
    # Arms Module
    Arms
        # # Methods
        # pull!, tick!, reset!,
        # # Arm Models
        # ArmBase,
        # # Stationary Arm Models
        # StationaryArmBase,
        # # Non Staionary Arm Models
        # NonstationaryArmBase

end # module


text:
using Parameters

mutable struct Region
	candidates::Int
	stars::Int
end

mutable struct Cell
	row::Int
	col::Int
	group::Int
	is_star::Bool
	is_blocked::Bool
end

@with_kw mutable struct Node
	data::Array{Cell, 2}
	cell::Union{Cell, Nothing} = nothing
	parent::Union{Node, Nothing} = nothing
	children::Vector{Node} = []
end

function insert!(n::Node, cells::Array{Cell}, cell::Cell)
	m = Node(data=deepcopy(cells), cell=cell, parent=n)
	push!(n.children, m)
end

function delete_children!(n::Node)
	for child in n.children
		delete_children!(child)
		deleteat!(n.children, findall(x->x == child, n.children))
	end
end


text:
function iterator end
macro iterator(ProxyType)
  ProxyType = esc(ProxyType)
  quote
    Base.iterate(p::$ProxyType) = Base.iterate(ProxyInterfaces.iterator(p))
    Base.iterate(p::$ProxyType, state) = Base.iterate(ProxyInterfaces.iterator(p), state)
    Base.IteratorSize(T::Type{<:$ProxyType}) = Base.IteratorSize(ProxyInterfaces.iterator(T))
    Base.IteratorEltype(T::Type{<:$ProxyType}) = Base.IteratorEltype(ProxyInterfaces.iterator(T))
    Base.eltype(T::Type{<:$ProxyType}) = Base.eltype(ProxyInterfaces.iterator(T))
    Base.length(p::$ProxyType) = Base.length(ProxyInterfaces.iterator(p))
    Base.size(p::$ProxyType) = Base.size(ProxyInterfaces.iterator(p))
    Base.size(p::$ProxyType, d) = Base.size(ProxyInterfaces.iterator(p), d)
    Base.axes(p::$ProxyType) = Base.axes(ProxyInterfaces.iterator(p))  # analog to Base.Generator
    Base.ndims(p::$ProxyType) = Base.ndims(ProxyInterfaces.iterator(p))
    # we include foreach to iterator, as it directly corresponds to a for loop
    Base.foreach(f, p::$ProxyType) = Base.foreach(f, ProxyInterfaces.iterator(p))
    Base.map(f, p::$ProxyType) = $ProxyType(Base.map(f, ProxyInterfaces.iterator(p)))
  end
end


text:
using Distributions

include("pdgsbpUtils.jl")
include("genericPDGSBP.jl")
include("pddpUtils.jl")
include("genericPDDP.jl")

x = toy_data([100, 50, 100], [-30 -20 -10;-20 0. 20;-30 10 20], ones(3,3), 1/3 * ones(3,3))

(maxiter, burnin, seed) = (150000, 20000, 1)
precompile(genericPDGSBP, (maxiter, burnin, seed))
precompile(genericPDDP, (maxiter, burnin, seed))

# run the PDGSBP sampler
@time predictive, P, clusters, times = genericPDGSBP(x, maxiter, burnin, seed, "/example_res")

max_stars_repo_name: cmerkatas/PDGSBP

max_issues_repo_name: cmerkatas/PDGSBP

max_forks_repo_name: cmerkatas/PDGSBP


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

using Markdown
using InteractiveUtils

# ╔═╡ 1971addc-cd61-4449-831e-8e3a4b449777
using PlutoTest

# ╔═╡ b90cb740-b02c-11eb-3bda-67e10d023147
md"""
# [Day 10: Elves Look, Elves Say](https://adventofcode.com/2015/day/10)

Today, the Elves are playing a game called [look-and-say](https://en.wikipedia.org/wiki/Look-and-say_sequence). They take turns making sequences by reading aloud the previous sequence and using that reading as the next sequence. For example, `211` is read as "one two, two ones", which becomes `1221` (`1` `2`, `2` `1`s).

Look-and-say sequences are generated iteratively, using the previous value as input for the next step. For each step, take the previous value, and replace each run of digits (like `111`) with the number of digits (`3`) followed by the digit itself (`1`).

For example:

-   `1` becomes `11` (`1` copy of digit `1`).
-   `11` becomes `21` (`2` copies of digit `1`).
-   `21` becomes `1211` (one `2` followed by one `1`).
-   `1211` becomes `111221` (one `1`, one `2`, and two `1`s).
-   `111221` becomes `312211` (three `1`s, two `2`s, and one `1`).

Starting with the digits in your puzzle input, apply this process 40 times. **What is the length of the result?**
"""

# ╔═╡ 0353fa91-ea22-4446-a83d-bb43cb52048d
puzzle_input = "3113322113"

# ╔═╡ 499c0130-654b-4d21-8f2f-1634b7e2b808
function groupsplit(s)
	groups = []
	last = nothing
	last_index = 1
	for idx in 1:length(s)
		if s[idx] != last
			if last != nothing
				push!(groups, s[last_index:idx-1])
			end
			last = s[idx]
			last_index = idx
		end
	end
	if last != nothing
		push!(groups, s[last_index:end])
	end

	groups
end

# ╔═╡ 32aabdc4-e704-4f7c-b3bd-89b76a9be3ee
lookandsay(s) = map(g -> string(length(g)) * g[1], groupsplit(s)) |> join

# ╔═╡ 3c013416-13bf-4113-889b-b37b3f6bbd15
lookandsay("11")

# ╔═╡ cc43c879-6720-4b47-906a-d37486fa5ea9
@assert lookandsay("1") == "11"

# ╔═╡ 547b7985-ed5b-4b25-93cd-21e05c54a563
@assert lookandsay("11") == "21"

# ╔═╡ 98f18475-0269-40e4-9cb8-85cb939a2b7d
@assert lookandsay("21") == "1211"

# ╔═╡ 2340661d-4f63-45ce-b4c9-a1cd314443ff
@assert lookandsay("1211") == "111221"

# ╔═╡ 0287d0bf-0870-4430-bb52-c7f6771b97dd
@assert lookandsay("111221") == "312211"

# ╔═╡ 43925c36-f5db-4895-9c6f-59c583af99b1
function applylookandsay(s, count)
	for i in 1:count
		s = lookandsay(s)
	end
	s
end

# ╔═╡ 8daceef8-005e-49dd-b9b5-5aed9e1bf83b
@assert applylookandsay("1", 5) == "312211"

# ╔═╡ 82717025-5412-4f8d-b0be-40b3bbac453a
part1 = length(applylookandsay(puzzle_input, 40))

# ╔═╡ 6bf5cb81-cb0b-45c0-967c-c3b6e8e46027
@test part1 == 329356

# ╔═╡ f113e0d2-fee5-40d8-96bf-0095d6d07eeb
md"Your puzzle answer was `329356`."

# ╔═╡ f264f5a9-5f6b-4604-80db-64a07f5f38c3
md"""
# Part Two

Neat, right? You might also enjoy hearing [John Conway talking about this sequence](https://www.youtube.com/watch?v=ea7lJkEhytA) (that's Conway of Conway's Game of Life fame).

Now, starting again with the digits in your puzzle input, apply this process 50 times. **What is the length of the new result?**
"""

# ╔═╡ 052995b1-fb2d-4ec5-92f4-63d707aa7e7d
part2 = length(applylookandsay(puzzle_input, 50))

# ╔═╡ b3d9ca1e-6757-4701-8658-accd25c033bf
@test part2 == 4666278

# ╔═╡ 768ff05c-909d-4738-bce8-949c8c96541f
md"Your puzzle answer was `4666278`."

# ╔═╡ Cell order:
# ╠═1971addc-cd61-4449-831e-8e3a4b449777
# ╟─b90cb740-b02c-11eb-3bda-67e10d023147
# ╟─0353fa91-ea22-4446-a83d-bb43cb52048d
# ╠═499c0130-654b-4d21-8f2f-1634b7e2b808
# ╠═32aabdc4-e704-4f7c-b3bd-89b76a9be3ee
# ╠═3c013416-13bf-4113-889b-b37b3f6bbd15
# ╠═cc43c879-6720-4b47-906a-d37486fa5ea9
# ╠═547b7985-ed5b-4b25-93cd-21e05c54a563
# ╠═98f18475-0269-40e4-9cb8-85cb939a2b7d
# ╠═2340661d-4f63-45ce-b4c9-a1cd314443ff
# ╠═0287d0bf-0870-4430-bb52-c7f6771b97dd
# ╠═43925c36-f5db-4895-9c6f-59c583af99b1
# ╠═8daceef8-005e-49dd-b9b5-5aed9e1bf83b
# ╠═82717025-5412-4f8d-b0be-40b3bbac453a
# ╠═6bf5cb81-cb0b-45c0-967c-c3b6e8e46027
# ╟─f113e0d2-fee5-40d8-96bf-0095d6d07eeb
# ╟─f264f5a9-5f6b-4604-80db-64a07f5f38c3
# ╠═052995b1-fb2d-4ec5-92f4-63d707aa7e7d
# ╠═b3d9ca1e-6757-4701-8658-accd25c033bf
# ╟─768ff05c-909d-4738-bce8-949c8c96541f


text:
#
#
#   These methods are for testing purposes only.
#   They are not optimized at all.
#
#

"""
Calculates `B_M B_M-1 ... B_1` and stabilizes the
matrix products by intermediate matrix decompositions.
Assumes that the input `Bs` are ordered as `[B_1, B_2, ..., B_M]`.
"""
function calc_Bchain_svd(Bs; svdalg = svd)
    M = length(Bs)
    svs = zeros(Float64, M,size(Bs[1],1))
    F = svdalg(Bs[1])
    U, S, Vt = F.U, F.S, F.Vt
    svs[1,:] = log.(S)
    svc = 2

    for k in 2:M
        # multiply B from the left to USVt
        # and update USVt
        U = Bs[k] * U
        U *= Diagonal(S)
        M = svdalg(U)
        U, S, Vtnew = M.U, M.S, M.Vt
        Vt = Vtnew * Vt

        # keep singular values
        svs[svc,:] = log.(S)
        svc += 1
    end

    return SVD(U, S, Vt), svs
end

function calc_Bchain_svd_faster(Bs; svdalg = svd)
    M = length(Bs)
    flv = 4
    N = 100
    G = eltype(eltype(Bs))

    start = 1
    stop = M

    U = Matrix{G}(I, flv*N, flv*N)
    S = ones(Float64, flv*N)
    Vt = Matrix{G}(I, flv*N, flv*N)

    tmp = copy(U)
    tmp2 = copy(Vt)

    svs = zeros(length(start:stop),flv*N)
    svc = 1
    for k in start:stop
        if mod(k,1) == 0 || k == stop # always decompose in the end
            mul!(tmp, Bs[k], U) # U = Bs[k] * U
            rmul!(tmp, Diagonal(S))
            M = svdalg(tmp)
            U, S, Vtnew = M.U, M.S, M.Vt
            mul!(tmp2, Vtnew, Vt)
            Vt .= tmp2

            svs[svc,:] = log.(S)
            svc += 1
        else
            U = Bs[k] * U
        end
    end
    return SVD(U, S, Vt), svs
end

"""
Calculates `B_M B_M-1 ... B_1` and stabilizes the
matrix products by intermediate matrix decompositions.
Assumes that the input `Bs` are ordered as `[B_1, B_2, ..., B_M]`.
"""
function calc_Bchain_qr(Bs)
    M = length(Bs)
    svs = zeros(Float64, M,size(Bs[1],1))
    U, D, T = udt(Bs[1])
    svs[1,:] = log.(D)
    svc = 2

    for k in 2:M
        # multiply B from the left to UDT
        # and update UDT
        U = Bs[k] * U
        U *= Diagonal(D)
        U, D, Tnew = udt(U)
        T = Tnew * T

        # keep singular values
        svs[svc,:] = log.(D)
        svc += 1
    end

    return UDT(U, D, T), svs
end

# """
# Calculate fake "G(τ, 0)", i.e. [B^-N1 + B^N2]^-1
# """
# function calc_tdgf_qr(B, N1, N2; inv_sum_method = inv_sum)
#   if N1 != 0
#     udtl = calc_Bchain_qr(inv(B), N1)[1]
#   else
#     udtl = udt!(Matrix(I*one(eltype(B)), size(B)))
#   end
#
#   if N2 != 0
#     udtr = calc_Bchain_qr(B, N2)[1]
#   else
#     udtr = udt!(Matrix(I*one(eltype(B)), size(B)))
#   end
#
#   inv_sum_method(udtl, udtr)
# end
#
#
# """
# Calculate fake "G(τ, 0)", i.e. [B^-N1 + B^N2]^-1
# """
# function calc_tdgf_svd(B, N1, N2; inv_sum_method = inv_sum, svdalg_inv = svd!, svdalg_chain = svd)
#   # svdalg_nonmutating = Symbol(replace(string(nameof(svdalg)), "!" => ""))
#   if N1 != 0
#     svdl = calc_Bchain_svd(inv(B), N1; svdalg = svdalg_chain)[1]
#   else
#     svdl = svdalg(Matrix(I*one(eltype(B)), size(B)))
#   end
#
#   if N2 != 0
#     svdr = calc_Bchain_svd(B, N2; svdalg = svdalg_chain)[1]
#   else
#     svdr = svdalg(Matrix(I*one(eltype(B)), size(B)))
#   end
#
#   inv_sum_method(svdl, svdr; svdalg = svdalg_inv)
# end


text:
module TestMacrosOrderBy

import Test: @testset, @test
import Volcanito: @order_by, materialize
import DataFrames: DataFrame

@testset "@order_by" begin
    df = DataFrame(
        a = [1, 2, 3, 4],
        b = [10, 20, 20, 10],
    )

    plan = @order_by(df, b)
    res = materialize(plan)
    @test res == DataFrame(
        a = [1, 4, 2, 3],
        b = [10, 10, 20, 20],
    )
end

end


text:
function parse_ref(ref_file::String, chroms::UnitRange)
    df = CSV.File(ref_file; types=Dict(:A1=>Char,:A2=>Char)) |> DataFrame
    df.A1 = tochar.(df.A1)
    df.A2 = tochar.(df.A2)
    @assert df.CHR isa Vector{Int}
    @assert df.BP isa Vector{Int}
    @assert df.MAF isa Vector{T} where T<:Real
    filter!(row->row.CHR in chroms, df)
    return df
end
parse_ref(ref_file::String, chrom::Integer) =
    parse_ref(ref_file, chrom:chrom)
function tochar(x)::Union{Char,Missing}
    if x isa String
        if length(x) == 1
            return first.(x)
        else
            return missing
        end
    else
        x
    end
end

function parse_bim(bim_file::String, chroms::UnitRange)
    header = [:CHR, :SNP, :POS, :BP, :A1, :A2]
    df = CSV.File(bim_file*".bim"; header=header, types=Dict(:A1=>Char,:A2=>Char)) |> DataFrame
    df.A1 = tochar.(df.A1)
    df.A2 = tochar.(df.A2)
    @assert df.CHR isa Vector{Int}
    filter!(row->row.CHR in chroms, df)
    return df
end
parse_bim(bim_file::String, chrom::Integer) =
    parse_bim(bim_file, chrom:chrom)

nuc_map(char::Char) = nuc_map(Val(char))
nuc_map(::Val{'A'}) = 'T'
nuc_map(::Val{'T'}) = 'A'
nuc_map(::Val{'C'}) = 'G'
nuc_map(::Val{'G'}) = 'C'
function permute_snps(df)
    unique(vcat(
        df[:,[:SNP,:A1,:A2]],
        DataFrame(SNP=df.SNP,
                  A1=df.A2,
                  A2=df.A1),
        DataFrame(SNP=df.SNP,
                  A1=nuc_map.(first.(df.A1)),
                  A2=nuc_map.(first.(df.A2))),
        DataFrame(SNP=df.SNP,
                  A1=nuc_map.(first.(df.A2)),
                  A2=nuc_map.(first.(df.A1)))
    ))
end
function join_snps(ref_df, vld_df, sst_df; verbose=false)
    # TODO: Be more efficient, don't allocate all this memory
    vld_snps = vld_df[:,[:SNP,:A1,:A2]]
    ref_snps = permute_snps(ref_df)
    sst_snps = permute_snps(sst_df)
    snps = innerjoin(vld_snps, ref_snps, sst_snps, on=[:SNP,:A1,:A2], makeunique=true)
    verbose && @info "$(nrow(snps)) common SNPs"
    return snps
end
norm_ppf(x) = quantile(Normal(), x)
function parse_sumstats(ref_df, vld_df, sst_file, n_subj; verbose=false, missingstring="")
    sst_df = CSV.File(sst_file; missingstring=missingstring, types=Dict(:A1=>Char,:A2=>Char)) |> DataFrame
    sst_df.A1 = tochar.(sst_df.A1)
    sst_df.A2 = tochar.(sst_df.A2)
    nucs = Set(['A','C','T','G'])
    filter!(row->(row.A1 in nucs) && (row.A2 in nucs), sst_df)
    filter!(row->!(row.P isa Missing) && !(row.BETA isa Missing), sst_df)
    sst_df.P = convert(Vector{Float64}, sst_df.P)
    sst_df.BETA = convert(Vector{Float64}, sst_df.BETA)
    snps = join_snps(ref_df, vld_df, sst_df; verbose=verbose)
    sort!(snps, [:SNP, :A1, :A2])

    n_sqrt = sqrt(n_subj)
    sst_eff = Dict{String,Float64}()
    for row in Tables.namedtupleiterator(sst_df)
        if hassnp(snps, (row.SNP,row.A1,row.A2)) ||
           hassnp(snps, (row.SNP,nuc_map.(row.A1),nuc_map.(row.A2)))
            effect_sign = 1
        elseif hassnp(snps, (row.SNP,row.A2,row.A1)) ||
               hassnp(snps, (row.SNP,nuc_map.(row.A2),nuc_map.(row.A1)))
            effect_sign = -1
        else
            continue
        end
        if hasproperty(row, :BETA)
            beta = row.BETA
        elseif hasproperty(row, :OR)
            beta = log(row.OR)
        end
        p = max(row.P, 1e-323)
        beta_std = effect_sign*sign(beta)*abs(norm_ppf(p/2))/n_sqrt
        sst_eff[row.SNP] = beta_std
    end
    _sst_df = DataFrame(SNP=String[],CHR=Int[],BP=Int[],BETA=Float64[],A1=Char[],A2=Char[],MAF=Float64[],FLP=Int[])
    for (idx,row) in enumerate(Tables.namedtupleiterator(ref_df))
        haskey(sst_eff, row.SNP) || continue

        SNP = row.SNP
        CHR = row.CHR
        BP = row.BP
        BETA = sst_eff[row.SNP]
        A1,A2 = row.A1,row.A2
        if hassnp(snps, (SNP,A1,A2))
            MAF = row.MAF
            FLP = 1
        elseif hassnp(snps, (SNP,A2,A1))
            A1, A2 = A2, A1
            MAF = 1-row.MAF
            FLP = -1
        elseif hassnp(snps, (SNP,nuc_map(A1),nuc_map(A2)))
            A1, A2 = nuc_map(A1), nuc_map(A2)
            MAF = row.MAF
            FLP = 1
        elseif hassnp(snps, (SNP,nuc_map(A2),nuc_map(A1)))
            A1, A2 = nuc_map(A2), nuc_map(A1)
            MAF = 1-row.MAF
            FLP = -1
        else
            verbose && @warn "(Chromosome $CHR) Didn't find ($SNP,$A1,$A2) in snps"
            # FIXME: Skip?
        end
        push!(_sst_df, (SNP=SNP,CHR=CHR,BP=BP,BETA=BETA,A1=A1,A2=A2,MAF=MAF,FLP=FLP))
    end
    return _sst_df
end

function findsnp(snps, (snp,a1,a2))
    SNP_range = binary_range_search(snps, snp, :SNP)
    SNP_range === nothing && return nothing
    SNP_L, SNP_R = SNP_range
    SNP_sub = snps[SNP_L:SNP_R,:]

    A1_range = binary_range_search(SNP_sub, a1, :A1)
    A1_range === nothing && return nothing
    A1_L, A1_R = A1_range
    A1_sub = SNP_sub[A1_L:A1_R,:]

    A2_range = binary_range_search(A1_sub, a2, :A2)
    A2_range === nothing && return nothing
    A2_L, A2_R = A2_range
    @assert A2_L == A2_R
    return SNP_L + (A1_L-1) + (A2_L-1)
end
hassnp(snps, row) = findsnp(snps, row) !== nothing
# TODO: Allow warm-restarts
function binary_range_search(snps, x, col)
    _snps = snps[!,col]
    L = 1
    R = nrow(snps)
    while true
        (L > R) && return nothing
        M = floor(Int, (L+R)/2)
        _x = _snps[M]
        if _x == x
            L,R = M,M
            snps_rows = nrow(snps)
            while (L > 1) && (_snps[L - 1] == x)
                L -= 1
            end
            while R < (snps_rows) && (_snps[R + 1] == x)
                R += 1
            end
            return L,R
        elseif _x < x
            L = M+1
        elseif _x > x
            R = M-1
        end
    end
end

function parse_ldblk(ldblk_dir, sst_df, chrom)
    chr_name = ldblk_dir * "/ldblk_1kg_chr" * string(chrom) * ".hdf5"
    hdf_chr = h5open(chr_name, "r")
    n_blk = length(hdf_chr)
    ld_blk = [read(hdf_chr["blk_"*string(blk)]["ldblk"]) for blk in 1:n_blk]

    snp_blk = Vector{String}[]
    for blk in 1:n_blk
        push!(snp_blk, read(hdf_chr["blk_"*string(blk)]["snplist"]))
    end

    blk_size = Int[]
    _ld_blk = Matrix{Float64}[]
    mm = 1
    for blk in 1:n_blk
        idx = [(ii,findfirst(s->s==snp, sst_df.SNP)) for (ii, snp) in enumerate(snp_blk[blk]) if snp in sst_df.SNP]
        push!(blk_size, length(idx))
        if !isempty(idx)
            idx_blk = mm:(mm+length(snp_blk[blk])-1)
            flip = [sst_df.FLP[jj] for jj in last.(idx)]
            flipM = flip' .* flip
            _blk = Matrix{Float64}(undef, length(idx), length(idx))
            P = collect(Iterators.product(first.(idx),first.(idx)))
            for icol in 1:size(P,2)
                for irow in 1:size(P,1)
                    row,col = P[irow,icol]
                    _blk[irow,icol] = ld_blk[blk][row,col] * flipM[irow,icol]
                end
            end
            push!(_ld_blk, _blk)
            mm += length(snp_blk[blk])
        else
            push!(_ld_blk,  Matrix{Float64}(undef, 0, 0))
        end
    end

    return _ld_blk, blk_size
end


text:
@muladd function perform_step!(integrator,
                               cache::Union{ISSEMConstantCache,
                                            ISSEulerHeunConstantCache},
                                            f=integrator.f)
  @unpack t,dt,uprev,u,p = integrator
  @unpack uf = cache
  nlsolve! = cache.nlsolve; nlcache = nlsolve!.cache
  alg = unwrap_alg(integrator, true)
  theta = alg.theta
  alg.symplectic ? a = dt/2 : a = dt

  # TODO: Stochastic extrapolants?
  u = uprev

  repeat_step = false
  if nlsolve! isa NLNewton
    uf.t = t
    J, nlcache.W = calc_W!(integrator, cache, dt*theta, repeat_step)
  end

  L = integrator.g(uprev,p,t)
  ftmp = integrator.f(uprev,p,t)

  if alg.symplectic
    z = zero(u) # constant extrapolation, justified by ODE IM
  else
    z = dt*ftmp # linear extrapolation
  end
  nlcache.z = z

  nlcache.c = a
  if alg.symplectic
    # u = uprev + z then  u = (uprev+u)/2 = (uprev+uprev+z)/2 = uprev + z/2
    #u = uprev + z/2
    tmp = uprev
  else
    tmp = uprev + dt*(1-theta)*ftmp
  end
  nlcache.tmp = tmp

  (z, η, iter, fail_convergence) = nlsolve!(integrator); fail_convergence && return nothing

  if alg.symplectic
    u = tmp + z
  else
    u = tmp + theta*z
  end

  gtmp = L.*integrator.W.dW

  if typeof(cache) <: ISSEulerHeunConstantCache
    utilde = u + gtmp
    gtmp = ((integrator.g(utilde,p,t) + L)/2)*integrator.W.dW
  end

  u += gtmp

  nlcache.ηold = η
  nlcache.nl_iters = iter

  if integrator.opts.adaptive

    Ed = dt*(J*ftmp)/2

    if typeof(cache) <: SplitStepEulerConstantCache
      K = @. uprev + dt * ftmp
      utilde =  @. K + integrator.sqdt * L
      ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
      En = ggprime .* (integrator.W.dW.^2 .- dt)./2
    elseif typeof(cache) <: ISSEulerHeunConstantCache
      utilde = uprev + L*integrator.sqdt
      ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
      En = ggprime.*(integrator.W.dW.^2)./2
    end

    resids = calculate_residuals(Ed, En, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(resids)
  end

  integrator.u = u
end

@muladd function perform_step!(integrator, cache::Union{ISSEMCache,
                                                        ISSEulerHeunCache},
                               f=integrator.f)
  @unpack t,dt,uprev,u,p = integrator
  @unpack uf,du1,dz,z,k,J,W,jac_config,gtmp,gtmp2,tmp,tmp,dW_cache = cache
  nlsolve! = cache.nlsolve; nlcache = nlsolve!.cache
  alg = unwrap_alg(integrator, true)
  alg.symplectic ? a = dt/2 : a = dt
  dW = integrator.W.dW
  mass_matrix = integrator.f.mass_matrix
  theta = alg.theta

  repeat_step = false

  if integrator.success_iter > 0 && !integrator.u_modified && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
  elseif alg.extrapolant == :linear
    @. u = uprev + integrator.fsalfirst*dt
  else # :constant
    copyto!(u,uprev)
  end

  nlsolve! isa NLNewton && calc_W!(integrator, cache, dt, repeat_step)

  integrator.f(tmp,uprev,p,t)

  if alg.symplectic
    @. z = zero(eltype(u)) # Justified by ODE solvers, constrant extrapolation when IM
  else
    @. z = dt*tmp # linear extrapolation
  end

  if alg.symplectic
    #@. u = uprev + z/2
    @. tmp = uprev
  else
    #@. u = uprev + dt*(1-theta)*tmp + theta*z
    @. tmp = uprev + dt*(1-theta)*tmp
  end
  nlcache.c = a
  (z, η, iter, fail_convergence) = nlsolve!(integrator); fail_convergence && return nothing

  if alg.symplectic
    @. u = uprev + z
  else
    #@. u = uprev + dt*(1-theta)*tmp + theta*z
    @. u = tmp + theta*z
  end

  nlcache.ηold = η
  nlcache.nl_iters = iter

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

  # Handle noise computations

  integrator.g(gtmp,uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    @. gtmp2 = gtmp*dW
  else
    mul!(gtmp2,gtmp,dW)
  end

  if typeof(cache) <: ISSEulerHeunCache
    gtmp3 = cache.gtmp3
    @. z = uprev + gtmp2
    integrator.g(gtmp3,z,p,t)
    @. gtmp = (gtmp3 + gtmp)/2
    if is_diagonal_noise(integrator.sol.prob)
      @. gtmp2 = gtmp*dW
    else
      mul!(gtmp2,gtmp,dW)
    end
  end

  @. u += gtmp2

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

  if integrator.opts.adaptive

    if has_invW(f)
      # This means the Jacobian was never computed!
      f.jac(J,uprev,p,t)
    end

    mul!(vec(z),J,vec(tmp))
    @. k = dt*dt*z/2

    # k is Ed
    # dz is En

      if !is_diagonal_noise(integrator.sol.prob)
        g_sized = norm(gtmp,2)
      else
        g_sized = gtmp
      end

      if typeof(cache) <: ISSEMCache
        @. z = uprev + dt*tmp + integrator.sqdt * g_sized

        if !is_diagonal_noise(integrator.sol.prob)
          integrator.g(gtmp,z,p,t)
          g_sized2 = norm(gtmp,2)
          @. dW_cache = dW.^2 - dt
          diff_tmp = integrator.opts.internalnorm(dW_cache)
          En = (g_sized2-g_sized)/(2integrator.sqdt)*diff_tmp
          @. dz = En
        else
          integrator.g(gtmp2,z,p,t)
          g_sized2 = gtmp2
          @. dz = (g_sized2-g_sized)/(2integrator.sqdt)*(dW.^2 - dt)
        end

      elseif typeof(cache) <: ISSEulerHeunCache
        @. z = uprev + integrator.sqdt * g_sized

        if !is_diagonal_noise(integrator.sol.prob)
          integrator.g(gtmp,z,p,t)
          g_sized2 = norm(gtmp,2)
          @. dW_cache = dW.^2
          diff_tmp = integrator.opts.internalnorm(dW_cache)
          En = (g_sized2-g_sized)/(2integrator.sqdt)*diff_tmp
          @. dz = En
        else
          integrator.g(gtmp2,z,p,t)
          g_sized2 = gtmp2
          @. dz = (g_sized2-g_sized)/(2integrator.sqdt)*(dW.^2)
        end

    end

    calculate_residuals!(tmp, k, dz, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(tmp)

  end
end


text:
module Arrangement
    using IntervalTrees
    using NearestNeighbors
    using Triangle
	using SparseArrays
	using LinearAlgebra
	using Distributed
	using DataStructures
	using ..Lar;
	using LinearAlgebraicRepresentation

	LA = LinearAlgebra

    include("./minimal_cycles.jl")
	include("./spatial_indices.jl")
    # include("./dimension_travel.jl")
    include("./planar_arrangement.jl")
    include("./spatial_arrangement.jl")

    # export planar_arrangement, spatial_arrangement

end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c54da25c-0dbf-41d1-bd7f-657b115bf6bf
md"""
# Test Notebook for [PlutoVista](https://github.com/j-fu/PlutVista.jl)
"""

# ╔═╡ 50d7298b-40ac-429c-8e4a-7e4d980bbfab
develop=true

# ╔═╡ 21b1f4b6-2172-4292-8249-35a615db4f59
begin
    using Pkg
    Pkg.activate(mktempdir())
	Pkg.add("Revise");using Revise
    Pkg.add("PlutoUI")
    if develop
     	Pkg.develop("PlutoVista")
    else
	    Pkg.add("PlutoVista")
    end
    using PlutoUI
    using PlutoVista
end

# ╔═╡ 825f6068-0f02-44e4-b083-1be3eb4c764d
let
	X=0:0.1:10
	p=PlutoCanvasPlot(;xrange=X,yrange=-1:1,resolution=(600,200))
	polyline!(p,X,sin.(X))
	polyline!(p,X,cos.(X))
    p
end

# ╔═╡ 2a5510d3-7355-4095-ab43-ad6d75bba090
let
	p=PlutoCanvasPlot(resolution=(600,300))

	axis!(p; xtics=(0:2.5:10),ytics=(-1:0.5:1))
	X=0:0.1:10
	linecolor!(p,1,0,0)
   	polyline!(p,X,sin.(X))
	linecolor!(p,0,1,0)
	polyline!(p,X,cos.(X))
end

# ╔═╡ 405b81de-3fd9-48f7-81a7-dc8f06dfda9c
let
	p=PlutoCanvasPlot(resolution=(300,300),xrange=-1:1,yrange=-1:1)
	T=0:0.01:2π
	X=0.45.*sin.(T)
	Y=0.45.*cos.(T)
	fillcolor!(p,0,1,0)
	polygon!(p,X,Y)
end

# ╔═╡ 83c16003-8196-4f2c-a602-c298f27475fb
begin
	@bind clock_t Clock(interval=0.01)
end

# ╔═╡ fd7c959b-94d6-43ff-b535-49b0883874d5
let
	p=PlutoCanvasPlot(resolution=(300,300),xrange=-1:1,yrange=-1:1)
	time=0.005*clock_t
	dt=0.01
    T=time:dt:2π/2+time
	px=5+sin(time/10)
	py=6+cos(time/10)
	X=sin.(px .*T)
	Y=cos.(py .*T)
	polyline!(p,X,Y)
end

# ╔═╡ Cell order:
# ╟─c54da25c-0dbf-41d1-bd7f-657b115bf6bf
# ╠═50d7298b-40ac-429c-8e4a-7e4d980bbfab
# ╠═21b1f4b6-2172-4292-8249-35a615db4f59
# ╠═825f6068-0f02-44e4-b083-1be3eb4c764d
# ╠═2a5510d3-7355-4095-ab43-ad6d75bba090
# ╠═405b81de-3fd9-48f7-81a7-dc8f06dfda9c
# ╠═83c16003-8196-4f2c-a602-c298f27475fb
# ╠═fd7c959b-94d6-43ff-b535-49b0883874d5


text:
using Test

# @testset "Core" begin
#   include("core/Core.jl")
# end

# @testset "Doctrines" begin
#   include("doctrines/Doctrines.jl")
# end

# @testset "WiringDiagrams" begin
#   include("wiring_diagrams/WiringDiagrams.jl")
# end

# @testset "Graphics" begin
#   include("graphics/Graphics.jl")
# end

# @testset "Programs" begin
#   include("programs/Programs.jl")
# end

@testset "LinearAlgebra" begin
  include("linear_algebra/LinearAlgebra.jl")
end


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_free_state_parameters(obj::AbstractGradientConverter)
    return jcall(obj, "getFreeStateParameters", jint, ())
end

function get_parameters(obj::DSSTGradientConverter, arg0::FieldSpacecraftState, arg1::DSSTForceModel)
    return jcall(obj, "getParameters", Vector{Gradient}, (FieldSpacecraftState, DSSTForceModel), arg0, arg1)
end

function get_state(obj::DSSTGradientConverter, arg0::DSSTForceModel)
    return jcall(obj, "getState", FieldSpacecraftState, (DSSTForceModel,), arg0)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
export WindowFunction

"""
    WindowFunction(_f::Function, width)

A radial window function.

It is defined by the callback `_f` between x = 0 (center) and x = `width`. x=1
is the position of the nearest neighbor and the range x ∈ [0,1] is called a
unit.

# Fields

- `_f<:Function`
Callback defining the window function.

It will be called as `_f(x::Real)` and is expected to return a `Real` number.
`x` can be expected to go from 0 (center of window) to `width`. `x` values above
`width` will return 0 without calling `_f`.

By convention, `_f(0) == 1`. If orthogonality to the nearest neighbor is
desired, it should return `_f(1) == 0`.

- `width<:Real`

Value up to which the callback function will be evaluated.

`width` must be bigger than 0. A `width` of 1 means that the window function is
evaluated up to the nearest neighbor.

# Examples

```jldoctest
julia> rectangular_window = WindowFunction(x -> 1., 0.5)
WindowFunction{var"#1#2", Float64}(var"#1#2"(), 0.5)

julia> rectangular_window(.4)
1.0

julia> rectangular_window(.6)
0.0
```
"""
struct WindowFunction{F<:Function, T<:Real}
    _f::F
    width::T

    function WindowFunction(_f::F, width::T) where {F<:Function, T<:Real}
        @assert width > 0 "WindowFunction width must be bigger than 0"
        new{F,T}(_f, width)
    end
end

function(w::WindowFunction)(x)
    x < 0 && error("undefined for negative values")
    x > w.width && return 0.0
    w._f(x)
end


text:
# Time to first plot

using BenchmarkTools

display("""include(joinpath("..", "src", "CairoPlot.jl")""")
@time include(joinpath("..", "src", "CairoPlot.jl"))
display("")

display("using Plots; gr()")
@time using Plots; gr()
display("using GR")
@time using GR
display("using Winston")
@time using Winston
display("using .CairoPlot")
@time using .CairoPlot

title = "Perturbation\n of a sine wave"
x = range(0, 10π, length=30)
y = sin.(x/5) .+ 0.3randn(length(x))
xticks = 1:5:30
yticks = -1.5:0.5:1.5

display("")
display("p = Plots.plot(x, y, title=title)")
@time p = Plots.plot(x, y, title=title)
display("g = GR.plot(x, y, title=title)")
@time g = GR.plot(x, y, title=title)
display("w = Winston.plot(x, y, title=title)")
@time w = Winston.plot(x, y, title=title)
display("c = CairoPlot.crplot(x, y; xticks, yticks, title)")
@time c = CairoPlot.crplot(x, y; xticks, yticks, title)

display("")
display("display(g)")
@time display(g)
display("display(p)")
@time display(p)
display("display(w)")
@time display(w)
display("display(c)")
@time display(c)


text:
"""
     Module for the kernel functions, also create kernel matrices
     Mostly from the list http://crsouza.com/2010/03/17/kernel-functions-for-machine-learning-applications/
     Kernels are created by calling the constructor (available right now :
     RBFKernel,
     Arguments are kernel specific, see for each one
"""
#LaplaceKernel, SigmoidKernel, ARDKernel, PolynomialKernel, Matern3_2Kernel, Matern5_2Kernel
module KernelModule

using LinearAlgebra
using Distances
using SpecialFunctions
using GradDescent: Optimizer
include("hyperparameters/HyperParametersModule.jl")
using .HyperParametersModule:
    Bound,
        OpenBound,
        ClosedBound,
        NullBound,
    Interval,
        interval,
    HyperParameter,
    HyperParameters,
        getvalue,
        setvalue!,
        checkvalue,
        gettheta,
        checktheta,
        settheta!,
        lowerboundtheta,
        upperboundtheta,
        update!,
        setfixed!,
        setfree!,
        setparamoptimizer!

import Base: *, +, getindex, show
export Kernel, KernelSum, KernelProduct
export RBFKernel, SEKernel, LaplaceKernel, SigmoidKernel, PolynomialKernel, ARDKernel, MaternKernel, Matern3_2Kernel, Matern5_2Kernel
export kernelmatrix,kernelmatrix!,kerneldiagmatrix,kerneldiagmatrix!
export computeIndPointsJ
export apply_gradients_lengthscale!, apply_gradients_variance!, apply_gradients!
export kernelderivativematrix_K
export kernelderivativematrix,kernelderivativediagmatrix
export InnerProduct, SquaredEuclidean, Identity
export compute_hyperparameter_gradient
export compute,plotkernel
export getvalue,setvalue!,setfixed!,setfree!,setoptimizer!, getvarianceoptimizer
export getlengthscales, getvariance
export isARD,isIso
export HyperParameter,HyperParameters,Interval, OpenBound,  NullBound

# include("KernelSum.jl")
# include("KernelProduct.jl")
include("kernel.jl")
include("RBF.jl")
include("matern.jl")
include("kernelmatrix.jl")
include("kernelgradients.jl")

end #end of module


text:
using CompTime
using Documenter

DocMeta.setdocmeta!(CompTime, :DocTestSetup, :(using CompTime); recursive=true)

makedocs(;
    modules=[CompTime],
    authors="Owen Lynch <root@owenlynch.org> and contributors",
    repo="https://github.com/olynch/CompTime.jl/blob/{commit}{path}#{line}",
    sitename="CompTime.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://olynch.github.io/CompTime.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/olynch/CompTime.jl",
    devbranch="main",
)


text:
using HDF5

# f = h5open("SARA/NatCom2020/results/outer_kernel_benchmark_energetic.h5")
# f = h5open("SARA/NatCom2020/results/outer_kernel_benchmark_uniform.h5")
f = h5open("SARA/NatCom2020/results/outer_kernel_benchmark_energetic_no_input_noise_short.h5")

using Plots
pyplot()

r2 = read(f, "likelihood ratio")
nsamples = read(f, "sample sizes")

tl = read(f, "temperature lengthscales")
dl = read(f, "dwelltime lengthscales")
names = read(f, "kernel names")

using Statistics
mr2 = median(r2, dims = 1)
mr2 = reshape(mr2, size(r2)[2:end])

# choose best length scale
br2 = zeros(size(mr2, 1), size(mr2, 2))
for i in eachindex(names)
    j = argmax(@view mr2[end, i, :, :])
    println(j)
    println(tl[j[1]])
    println(dl[j[2]])
    br2[:, i] = mr2[:, i, j[1], j[2]]
end
plot!(xlabel = "number of samples", ylabel = "R^2_sigma")
plot!(nsamples, br2, label = permutedims(names), legend = :bottomright)
gui()


text:
_subscript(i::Integer) = i<0 ? error("$i is negative") : join('₀'+d for d in reverse(digits(i)))

function _parameters(x::RiskParameters{M}) where {M <: EpidemicModel}
  inds = _indices(x, cumulative=false)
  params1 = ["ϵ", "Ωs", "κ", "Ωt", "Ωl", "Ωr"]
  params2 = String[]
  for k = eachindex(inds)
    for l = 1:inds[k]
      push!(params2, params1[k] * _subscript(l))
    end
  end
  return params2
end

"""
Produces a data frame summarizing MCMC results.
"""
function summary(mcmc::MCMC{M}; burnin::Int64=0, thin::Int64=1, bychain::Bool=true, credibleinterval::Float64=0.95) where M <: EpidemicModel
  !(0.0 < credibleinterval < 1.0) && error("`credibleinterval` must be between 0 and 1")
  ci = (1-credibleinterval)/2
  params = _parameters(mcmc.markov_chains[1].risk_parameters[1])
  if bychain && length(mcmc.markov_chains) > 1
    df = DataFrame(parameter=String[],
                   chain=Int64[],
                   mean=Float64[],
                   var=Float64[],
                   CI=Tuple{Float64, Float64}[])
    for i = 1:length(mcmc.markov_chains)
      trace_array = convert(Array{Float64, 2}, mcmc.markov_chains[i].risk_parameters[1+burnin:thin:end])
      append!(df,
        DataFrame(parameter = params,
                  chain = fill(i, size(trace_array, 2)),
                  mean = mean(trace_array, dims=1)[:],
                  var = var(trace_array, dims=1)[:],
                  CI = [Tuple(quantile(trace_array[:,j], [ci, 1-ci])) for j = 1:size(trace_array, 2)]))
    end
  else
    trace_array = convert(Array{Float64, 2}, mcmc.markov_chains[1].risk_parameters[1+burnin:thin:end])
    for i = 2:length(mcmc.markov_chains)
      vcat!(trace_array, convert(Array{Float64, 2}, mcmc.markov_chains[i].risk_parameters[1+burnin:thin:end]))
    end
    df = DataFrame(parameter = params,
                   mean = mean(trace_array, dims=1)[:],
                   var = var(trace_array, dims=1)[:],
                   CI = [Tuple(quantile(trace_array[:,j], [ci, 1-ci])) for j = 1:size(trace_array, 2)])
  end
  return df
end


text:
#= src/LowRankApprox.jl
=#
__precompile__()
module LowRankApprox
using FillArrays
using LinearAlgebra, SparseArrays, Random

import Base: convert,
             eltype, size, getindex, setindex!, copy,
             isreal, real, imag, Nothing, copyto!,
             +, -, *, /, \, ^, conj, conj!, promote_rule, similar, fill!,
             Matrix
import LinearAlgebra: BlasFloat, BlasInt, checksquare, chkstride1, rank

using FFTW
import FFTW: plan_r2r!, R2HC, r2rFFTWPlan, FFTWPlan
using Nullables
import LinearAlgebra: mul!, ldiv!, transpose, transpose!, axpy!, ishermitian, issymmetric,
            lmul!, rmul!, adjoint, adjoint!
import SparseArrays: sparse
ind2sub(dims, ind) = CartesianIndices(dims)[ind]
qrfact!(A) = qr!(A)

import FillArrays: AbstractFill

export

  # LowRankApprox.jl
  LRAOptions,

  # cur.jl
  AbstractCURPackedU, CURPackedU, HermitianCURPackedU, SymmetricCURPackedU,
  AbstractCUR, CUR, HermitianCUR, SymmetricCUR,
  curfact, curfact!, cur, cur!,

  # id.jl
  IDPackedV, ID, idfact, idfact!, id, id!,

  # linop.jl
  AbstractLinearOperator, AbstractMatOrLinOp,
  LinearOperator, HermitianLinearOperator,

  # matrixlib.jl
  matrixlib,

  # permute.jl
  PermutationMatrix, RowPermutation, ColumnPermutation,

  # pheig.jl
  HermitianPartialEigen, pheigfact, pheig, pheigvals,

  # pqr.jl
  PartialQRFactors, PartialQR, pqrfact, pqrfact!, pqr, pqr!,

  # prange.jl
  prange, prange!,

  # psvd.jl
  PartialSVD, psvdfact, psvd, psvdvals,

  # sketch.jl
  SketchMatrix, RandomGaussian, RandomSubset, SRFT, SparseRandomGaussian,
  sketch, sketchfact,

  # snorm.jl
  snorm, snormdiff,

  # trapezoidal.jl
  Trapezoidal, LowerTrapezoidal, UpperTrapezoidal,

  # aca.jl
  LowRankMatrix

# common

mutable struct LRAOptions
  atol::Float64
  maxdet_niter::Int
  maxdet_tol::Float64
  nb::Int
  pheig_orthtol::Float64
  pqrfact_retval::String
  rank::Int
  rtol::Float64
  sketch::Symbol
  sketch_randn_niter::Int
  sketchfact_adap::Bool
  sketchfact_randn_samp::Function
  sketchfact_srft_samp::Function
  sketchfact_sub_samp::Function
  snorm_niter::Int
  verb::Bool
end

function LRAOptions(::Type{T}; args...) where T
  opts = LRAOptions(
    0,                  # atol
    -1,                 # maxdet_niter
    -1,                 # maxdet_tol
    32,                 # nb
    sqrt(eps(real(T))), # pheig_orthtol
    "qr",               # pqrfact_retval
    -1,                 # rank
    5*eps(real(T)),     # rtol
    :randn,             # sketch
    0,                  # sketch_randn_niter
    true,               # sketchfact_adap
    n -> n + 8,         # sketchfact_randn_samp
    n -> n + 8,         # sketchfact_srft_samp
    n -> 4*n + 8,       # sketchfact_sub_samp
    32,                 # snorm_niter
    true,               # verb
  )
  for (key, value) in args
    setfield!(opts, key, value)
  end
  opts
end
LRAOptions(; args...) = LRAOptions(Float64; args...)

function copy(opts::LRAOptions; args...)
  opts_ = LRAOptions()
  for field in fieldnames(typeof(opts))
    setfield!(opts_, field, getfield(opts, field))
  end
  for (key, value) in args
    setfield!(opts_, key, value)
  end
  opts_
end

function chkopts!(opts::LRAOptions)
  opts.atol >= 0 || throw(ArgumentError("atol"))
  opts.nb > 0 || throw(ArgumentError("nb"))
  opts.pheig_orthtol >= 0 || throw(ArgumentError("pheig_orthtol"))
  opts.rtol >= 0 || throw(ArgumentError("rtol"))
  opts.sketch in (:none, :randn, :sprn, :srft, :sub) ||
    throw(ArgumentError("sketch"))
  opts.pqrfact_retval = lowercase(opts.pqrfact_retval)
end
function chkopts!(opts::LRAOptions, A)
  chkopts!(opts)
  if typeof(A) <: AbstractLinOp && opts.sketch != :randn
    warn("invalid sketch method; using \"randn\"")
    opts.sketch = :randn
  end
end

chktrans(trans::Symbol) = trans in (:n, :c) || throw(ArgumentError("trans"))

# source files

include("lapack.jl")
include("linop.jl")
include("matrixlib.jl")
include("permute.jl")
include("snorm.jl")
include("trapezoidal.jl")
include("util.jl")

include("cur.jl")
include("id.jl")
include("pheig.jl")
include("pqr.jl")
include("prange.jl")
include("psvd.jl")
include("sketch.jl")

include("lowrankmatrix.jl")

end  # module


text:
using LinearAlgebra

function powermethod(A, x0::AbstractVector=randn(size(A, 2)) |> normalize!; niter::Int=1000)
    e,v = eigen(A)
    return e[1], v[:,1]
    for i = 1:niter
        x0 = A*x0
        x0 |> normalize!
    end
    x0'*A*x0, x0
end

function powermethod_back(A, e, v, de, dv)
    @show A
    de isa Nothing && (de = zero(e))
    dv isa Nothing && (dv = zero(v))
    nm = norm(A*v)
    dvAinv = dv'*inv(A)
    Δ = zero(A)
    de*v*v' + norm(A*v)*(dvAinv*v)*v*v' - dvAinv'*v'
end

@primitive powermethod(A),dy,y powermethod_back(A, y..., dy...)

xa = randn(4,4); xa+=xa'
x = Param(xa)
b = randn(4)
y = @diff powermethod(x)[2] |> sum
grad(y, x)
#e = @diff powermethod(x)[1]
#grad(e, x)

δ = 1e-5
xa_ = copy(xa)
xa_[1]+=δ/2
pos = powermethod(xa_)
xa_[1]-=δ
neg = powermethod(xa_)
(pos[1]-neg[1])/δ

using Test
@testset "power method" begin
    N = 4
    A = randn(N, N); A = A+A'
    eg, vg = powermethod(A)
    eigenres = eigen(A)
    i = argmax(eigenres.values .|> abs)
    @test eg ≈ eigenres.values[i]
    @test vg'*eigenres.vectors[:,i] |> abs ≈ 1
    @test AutoGrad.gradcheck(A->powermethod(A)[2].*b |> sum, A)
end

N = 4
A = randn(N, N); A = A+A'
eg, vg = powermethod(A)
eigenres = eigen(A)
i = argmax(eigenres.values .|> abs)
@test eg ≈ eigenres.values[i]
@test vg'*eigenres.vectors[:,i] |> abs ≈ 1
@test AutoGrad.gradcheck(A->powermethod(A)[2].*b |> sum, A)


text:
module OptimNewton

function optimize_newton(x, f, g, h;
                         grtol=1e-4,
                         gatol=1e-4,
                         prints=false)
    while true
        fx = f(x)
        gx = g(x)
        hx = h(x)
        #convergence
        abs(abs(gx) / fx) <= grtol && break
        abs(gx) <= gatol && break

        if prints
            println("|G| = $(abs(gx))")
        end

        #move to next p
        x -= gx / hx
    end
    return x
end

end # module