text:
###############################################################################
#                                                                             #
#  This file is part of the julia module for Multi Objective Optimization     #
#  (c) Copyright 2017 by Aritra Pal                                           #
#                                                                             #
#  Permission is granted for academic research use.  For other uses,          #
#  contact the authors for licensing options.                                 #
#                                                                             #
#  Use at your own risk. I make no guarantees about the correctness or        #
#  usefulness of this code.                                                   #
#                                                                             #
###############################################################################

using Modof, DataFrames, DataFramesMeta, Match, PyCall, PyPlot

@pyimport seaborn as sns
sns.set_style("white")
sns.set_context("paper")
sns.set(font_scale=2.0)

plot_location = "../plots/"

#####################################################################
# Collecting the Experimental Results in one place                  #
#####################################################################

files = readdir("../results/")
files = files[[contains(files[i], "Experimental_Results") for i in 1:length(files)]]
data = readtable("../results/$(files[1])")
for i in 2:length(files)
	data = vcat(data, readtable("../results/$(files[i])"))
end
data = data[names(data)[1:21]]
data = data[completecases(data),:]
unique!(data)
sort!(data)
[rm("../results/$(files[i])") for i in 1:length(files)]
writetable("../results/Experimental_Results.csv", data)

#####################################################################
# Computing the results for AUTO for a Single Threads               #
#####################################################################

instances = unique(data[[:Problem_Type, :Instance_Type, :Instance_ID, :Correlation, :Objectives, :Variables, :Constraints, :Sparsity]])
algorithms = unique(data[[:Algorithm, :Type_of_Operation, :MIP_Emphasis, :Type_of_Boost, :Cplex_Heur, :Total_Threads]])
instance_inds = Int64[]
algorithm_inds = Int64[]
for i in 1:nrow(instances)
    push!(instance_inds, i)
    if (1-instances[:Sparsity][i]) >= 0.42
        push!(algorithm_inds, 10)
    else
        if floor(abs(log10(instances[:Constraints][i]/instances[:Variables][i]))) > 1.0
            push!(algorithm_inds, 4)
        else
            push!(algorithm_inds, 14)
        end
    end
end
results = hcat(instances[instance_inds,:], algorithms[algorithm_inds, :])
data_ = zeros(nrow(results), 5)
count = 1
for i in 1:nrow(results)
	tmp = data[data[:Problem_Type] .== results[:Problem_Type][i],:]
    tmp = tmp[tmp[:Instance_Type] .== results[:Instance_Type][i],:]
    tmp = tmp[tmp[:Instance_ID] .== results[:Instance_ID][i],:]
    tmp = tmp[tmp[:Algorithm] .== results[:Algorithm][i],:]
    tmp = tmp[tmp[:Type_of_Operation] .== results[:Type_of_Operation][i],:]
    tmp = tmp[tmp[:MIP_Emphasis] .== results[:MIP_Emphasis][i],:]
    tmp = tmp[tmp[:Type_of_Boost] .== results[:Type_of_Boost][i],:]
    tmp = tmp[tmp[:Cplex_Heur] .== results[:Cplex_Heur][i],:]
    tmp = tmp[tmp[:Total_Threads] .== results[:Total_Threads][i],:]
	if nrow(tmp) == 1
		data_[count, 1] = tmp[:Master_Time][1]
        data_[count, 2] = tmp[:Slave_Time][1]
        data_[count, 3] = tmp[:IP_Solved][1]
        data_[count, 4] = tmp[:Num_Non_Dom_Pts][1]
        data_[count, 5] = tmp[:Num_Non_Dom_Pts_Found][1]
	end
	count += 1
end
results[:Algorithm] = "AUTO"
results[:Type_of_Operation] = "lex_min"
results[:MIP_Emphasis] = false
results[:Type_of_Boost] = "NB"
results[:Cplex_Heur] = false
results[:Parallelism] = false
results[:Thread_ID] = 1
results[:Master_Time] = data_[:, 1]
results[:Slave_Time] = data_[:, 2]
results[:IP_Solved] = data_[:, 3]
results[:Num_Non_Dom_Pts] = data_[:, 4]
results[:Num_Non_Dom_Pts_Found] = data_[:, 5]
results = results[names(results)[[1:13..., 15, 14, 16:21...]]]
writetable("../results/Experimental_Results_0.csv", results)

#####################################################################
# Recollecting the Experimental Results in one place                #
#####################################################################

files = readdir("../results/")
files = files[[contains(files[i], "Experimental_Results") for i in 1:length(files)]]
data = readtable("../results/$(files[1])")
for i in 2:length(files)
	data = vcat(data, readtable("../results/$(files[i])"))
end
data = data[names(data)[1:21]]
data = data[completecases(data),:]
unique!(data)
sort!(data)
[rm("../results/$(files[i])") for i in 1:length(files)]
writetable("../results/Experimental_Results.csv", data)

#####################################################################
# Collecting the Experimental Results in one place                  #
#####################################################################

instances = unique(data[[:Problem_Type, :Instance_Type, :Instance_ID, :Correlation, :Objectives, :Variables, :Constraints, :Sparsity]])
algorithms = unique(data[[:Algorithm, :Type_of_Operation, :MIP_Emphasis, :Type_of_Boost, :Cplex_Heur, :Total_Threads]])
instance_inds = Int64[]
algorithm_inds = Int64[]
data_ = zeros(nrow(instances)*nrow(algorithms), 8)
count = 1
for i in 1:nrow(instances)
    tmp = data[data[:Problem_Type] .== instances[:Problem_Type][i],:]
    tmp = tmp[tmp[:Instance_Type] .== instances[:Instance_Type][i],:]
    tmp = tmp[tmp[:Instance_ID] .== instances[:Instance_ID][i],:]
    for j in 1:nrow(algorithms)
    	k = algorithms[:Total_Threads][j]
        tmp2 = tmp[tmp[:Algorithm] .== algorithms[:Algorithm][j],:]
        tmp2 = tmp2[tmp2[:Type_of_Operation] .== algorithms[:Type_of_Operation][j],:]
        tmp2 = tmp2[tmp2[:MIP_Emphasis] .== algorithms[:MIP_Emphasis][j],:]
        tmp2 = tmp2[tmp2[:Type_of_Boost] .== algorithms[:Type_of_Boost][j],:]
        tmp2 = tmp2[tmp2[:Cplex_Heur] .== algorithms[:Cplex_Heur][j],:]
        tmp2 = tmp2[tmp2[:Total_Threads] .== algorithms[:Total_Threads][j],:]
        if nrow(tmp2) == k
            push!(instance_inds, i)
            push!(algorithm_inds, j)
            data_[count, 1] = maximum(tmp2[:Master_Time]) + maximum(tmp2[:Slave_Time])
            data_[count, 2:(2+k-1)] = [tmp2[:Slave_Time]...]
            data_[count, 2] += maximum(tmp2[:Master_Time])
            data_[count, 6] = sum(tmp2[:IP_Solved])
            data_[count, 7] = mean(tmp2[:Num_Non_Dom_Pts])
            data_[count, 8] = sum(tmp2[:Num_Non_Dom_Pts_Found]) - (k-1)
            count += 1
        end
    end
end
data_ = data_[1:length(instance_inds),:]
results = DataFrame()
results[:Problem_Type] = instances[:Problem_Type][instance_inds]
results[:Instance_Type] = instances[:Instance_Type][instance_inds]
results[:Instance_ID] = instances[:Instance_ID][instance_inds]
results[:Correlation] = instances[:Correlation][instance_inds]
results[:Objectives] = instances[:Objectives][instance_inds]
results[:Variables] = instances[:Variables][instance_inds]
results[:Constraints] = instances[:Constraints][instance_inds]
results[:Sparsity] = instances[:Sparsity][instance_inds]
results[:Algorithm] = algorithms[:Algorithm][algorithm_inds]
results[:Type_of_Operation] = algorithms[:Type_of_Operation][algorithm_inds]
results[:MIP_Emphasis] = algorithms[:MIP_Emphasis][algorithm_inds]
results[:Type_of_Boost] = algorithms[:Type_of_Boost][algorithm_inds]
results[:Cplex_Heur] = algorithms[:Cplex_Heur][algorithm_inds]
results[:Total_Threads] = algorithms[:Total_Threads][algorithm_inds]
results[:Total_Time] = round(data_[:,1],4)
results[:Thread_1_Slave_Time] = round(data_[:,2],4)
results[:Thread_2_Slave_Time] = round(data_[:,3],4)
results[:Thread_3_Slave_Time] = round(data_[:,4],4)
results[:Thread_4_Slave_Time] = round(data_[:,5],4)
results[:IP_Solved] = data_[:,6]
results[:Percent_of_Non_Dom_Pts_Found] = round((data_[:,8] ./ data_[:,7] ) * 100, 4)
sort!(results)
writetable("../results/UER1.csv", results)

#####################################################################
# Compressing Algorithms and Number of Threads                      #
# into a Single Column                                              #
#####################################################################

algorithm_ = String[]
for i in 1:nrow(results)
	push!(algorithm_, results[:Algorithm][i])
	if algorithm_[end] != "AUTO"
		algorithm_[end] = @match results[:Type_of_Operation][i] begin
			"lex_min" => "$(algorithm_[end])_LM"
			"augmented_operation" => "$(algorithm_[end])_AO"
			"dynamic" => "$(algorithm_[end])_D"
		end
		algorithm_[end] = @match results[:MIP_Emphasis][i] begin
			true => "$(algorithm_[end])2"
			_ => algorithm_[end]
		end
		algorithm_[end] = @match results[:Type_of_Boost][i] begin
			"NB" => "$(algorithm_[end])_NB"
			"B" => "$(algorithm_[end])_B"
			"O_SH" => "$(algorithm_[end])_OSH"
			"M_SH" => "$(algorithm_[end])_MSH"
		end
		algorithm_[end] = @match results[:Cplex_Heur][i] begin
			true => "$(algorithm_[end])_HON"
			_ => "$(algorithm_[end])_HOFF"
		end
	end
	algorithm_[end] = "$(algorithm_[end])_$(results[:Total_Threads][i])"
end
results[:Algorithm] = algorithm_
delete!(results, [:Type_of_Operation, :MIP_Emphasis, :Type_of_Boost, :Cplex_Heur, :Total_Threads])
sort!(results)
writetable("../results/UER2.csv", results)

#####################################################################
# Computing the Rank of Algorithms for a Single Thread              #
#####################################################################

minimum_time = Float64[]
ratio_to_minimum_time = Float64[]
rank = Int64[]
for i in 1:nrow(results)
	tmp = results[results[:Problem_Type] .== results[:Problem_Type][i],:]
    tmp = tmp[tmp[:Instance_Type] .== results[:Instance_Type][i],:]
    tmp = tmp[tmp[:Instance_ID] .== results[:Instance_ID][i],:]
	thread = parse(Int64, split(results[:Algorithm][i], "_")[end])
	tmp = tmp[ [ contains(tmp[:Algorithm][j], "_$(thread)") for j in 1:nrow(tmp) ],:]
	push!(minimum_time, minimum(tmp[:Total_Time]))
	push!(ratio_to_minimum_time, results[:Total_Time][i]/minimum_time[end])
	tmp2 = sort(tmp[:Total_Time])
	push!(rank, findin(tmp2, results[:Total_Time][i])[1])
end
results[:Minimum_Time] = minimum_time
results[:Ratio_to_Minimum_Time] = ratio_to_minimum_time
results[:Rank_of_the_Algorithm] = rank
results = results[names(results)[[1:9..., 18, 19, 10:17...]]]
sort!(results)
writetable("../results/UER3.csv", results)

#####################################################################
# Decomposing Results for Single Threads and Multiple Threads       #
#####################################################################

single_thread_results = results[[ contains(results[:Algorithm][i], "_1") for i in 1:nrow(results) ], :]
for i in 1:nrow(single_thread_results)
	single_thread_results[:Algorithm][i] = @match single_thread_results[:Algorithm][i] begin
		"AUTO_1" => "AUTO"
		"UDECM_LM_NB_HON_1" => "V1"
 		"BBM_LM_NB_HON_1" => "V2"
 		"EBBM_LM_NB_HON_1" => "V3"
 		"UDECM_AO_NB_HON_1" => "V4"
 		"UDECM_D_NB_HON_1" => "V5"
 		"EBBM_LM2_NB_HON_1" => "V6"
 		"EBBM_AO_NB_HON_1" => "V7"
 		"EBBM_D_NB_HON_1" => "V8"
 		"UDECM_D_B_HON_1" => "V9"
 		"UDECM_D_OSH_HON_1" => "V10"
 		"EBBM_LM2_B_HON_1" => "V11"
 		"EBBM_LM2_OSH_HON_1" => "V12"
 		"EBBM_LM2_B_HOFF_1" => "V13"
 		"EBBM_LM2_OSH_HOFF_1" => "V14"
 		"EBBM_D_B_HON_1" => "V15"
 		"EBBM_D_OSH_HON_1" => "V16"
		_ => single_thread_results[:Algorithm][i]
	end
end
sort!(single_thread_results)
writetable("../results/UER4.csv", single_thread_results)

multiple_thread_results = results[[ contains(results[:Algorithm][i], "AUTO") for i in 1:nrow(results) ], :]
multiple_thread_results[:Algorithm] = [ split(multiple_thread_results[:Algorithm][i], "_")[end] for i in 1:nrow(multiple_thread_results)]
sort!(multiple_thread_results)
writetable("../results/UER5.csv", multiple_thread_results)

#####################################################################
# Decomposing Results for Small and Large Scale Instances           #
#####################################################################

large_instance_minimum_time = 600.0

single_thread_results_of_small_instances = @where(single_thread_results, :Minimum_Time .< large_instance_minimum_time)
sort!(single_thread_results_of_small_instances)
writetable("../results/UER4_SI.csv", single_thread_results_of_small_instances)
multiple_thread_results_of_small_instances = @where(multiple_thread_results, :Minimum_Time .< large_instance_minimum_time)
sort!(multiple_thread_results_of_small_instances)
writetable("../results/UER5_SI.csv", multiple_thread_results_of_small_instances)

single_thread_results_of_large_instances = @where(single_thread_results, :Minimum_Time .>= large_instance_minimum_time)
sort!(single_thread_results_of_large_instances)
writetable("../results/UER4_LI.csv", single_thread_results_of_large_instances)
multiple_thread_results_of_large_instances = @where(multiple_thread_results, :Minimum_Time .>= large_instance_minimum_time)
sort!(multiple_thread_results_of_large_instances)
writetable("../results/UER5_LI.csv", multiple_thread_results_of_large_instances)

#####################################################################
# Generating the DataFrame for additional insights                  #
#####################################################################

data = zeros(17, 15)
algorithm = String[]
for i in 1:17
	if i <= 16
		push!(algorithm, "V$i")
		tmp = @where(single_thread_results, :Algorithm .== "V$i")
	else
		push!(algorithm, "AUTO")
		tmp = @where(single_thread_results, :Algorithm .== "AUTO")
	end
	data[i, 1] = mean(tmp[:Total_Time])
	data[i, 2] = mean(tmp[:Ratio_to_Minimum_Time])
	data[i, 3] = mean(tmp[:Rank_of_the_Algorithm])
	data[i, 4] = mean(tmp[:IP_Solved])
	data[i, 5] = mean(tmp[:Percent_of_Non_Dom_Pts_Found])
	if i <= 16
		tmp = @where(single_thread_results_of_small_instances, :Algorithm .== "V$i")
	else
		tmp = @where(single_thread_results_of_small_instances, :Algorithm .== "AUTO")
	end
	data[i, 6] = mean(tmp[:Total_Time])
	data[i, 7] = mean(tmp[:Ratio_to_Minimum_Time])
	data[i, 8] = mean(tmp[:Rank_of_the_Algorithm])
	data[i, 9] = mean(tmp[:IP_Solved])
	data[i, 10] = mean(tmp[:Percent_of_Non_Dom_Pts_Found])
	if i <= 16
		tmp = @where(single_thread_results_of_large_instances, :Algorithm .== "V$i")
	else
		tmp = @where(single_thread_results_of_large_instances, :Algorithm .== "AUTO")
	end
	data[i, 11] = mean(tmp[:Total_Time])
	data[i, 12] = mean(tmp[:Ratio_to_Minimum_Time])
	data[i, 13] = mean(tmp[:Rank_of_the_Algorithm])
	data[i, 14] = mean(tmp[:IP_Solved])
	data[i, 15] = mean(tmp[:Percent_of_Non_Dom_Pts_Found])
end
results = DataFrame()
results[:Algorithm] = algorithm
results[:All_Instances_Total_Time] = data[:, 1]
results[:All_Instances_Ratio] = data[:, 2]
results[:All_Instances_Rank] = data[:, 3]
results[:All_Instances_IP_Solved] = data[:, 4]
results[:All_Instances_Percent] = data[:, 5]
results[:Small_Instances_Total_Time] = data[:, 6]
results[:Small_Instances_Ratio] = data[:, 7]
results[:Small_Instances_Rank] = data[:, 8]
results[:Small_Instances_IP_Solved] = data[:, 9]
results[:Small_Instances_Percent] = data[:, 10]
results[:Large_Instances_Total_Time] = data[:, 11]
results[:Large_Instances_Ratio] = data[:, 12]
results[:Large_Instances_Rank] = data[:, 13]
results[:Large_Instances_IP_Solved] = data[:, 14]
results[:Large_Instances_Percent] = data[:, 15]
writetable("../results/UER6.csv", results)

#####################################################################
# Generating the Plots                                              #
#####################################################################

#####################################################################
## Plots for the 3 classes of Instances                            ##
#####################################################################

#####################################################################
### Plots for All Instances                                       ###
#####################################################################

plt_performance_profiles(single_thread_results, [:Problem_Type, :Instance_Type, :Instance_ID], ["AUTO", "V1", "V2", "V3"], "../plots/1_Thread_All_Instances_Performance_Profiles.eps")

single_thread_results = single_thread_results[[single_thread_results[:Algorithm][i] in ["AUTO", "V4", "V5", "V6", "V7", "V8", "V10", "V15", "V16"] for i in 1:nrow(single_thread_results)],:]

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results[:Algorithm], y=single_thread_results[:Total_Time])
ax = gca()
ax[:set_yscale]("log")
xlabel("Algorithms")
ylabel("Computing Time (secs)")
title("All Instances")
f[:savefig]("../plots/1_Thread_All_Instances_Computing_Time.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results[:Algorithm], y=single_thread_results[:Ratio_to_Minimum_Time])
ax = gca()
ax[:set_yscale]("log")
xlabel("Algorithms")
ylabel("Ratio to Minimum Time")
title("All Instances")
f[:savefig]("../plots/1_Thread_All_Instances_Ratio_to_Minimum_Time.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results[:Algorithm], y=single_thread_results[:Rank_of_the_Algorithm])
xlabel("Algorithms")
ylabel("Rank of the Algorithm")
title("All Instances")
f[:savefig]("../plots/1_Thread_All_Instances_Rank_of_the_Algorithm.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results[:Algorithm], y=single_thread_results[:IP_Solved])
xlabel("Algorithms")
ylabel("IPs Solved")
title("All Instances")
f[:savefig]("../plots/1_Thread_All_Instances_IP_Solved.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results[:Algorithm], y=single_thread_results[:Percent_of_Non_Dom_Pts_Found])
xlabel("Algorithms")
ylabel("Percentage of Nondominated Points Found")
title("All Instances")
f[:savefig]("../plots/1_Thread_All_Instances_Percent_of_Non_Dom_Pts_Found.eps")

#####################################################################
### Plots for Small Instances                                     ###
#####################################################################

plt_performance_profiles(single_thread_results_of_small_instances, [:Problem_Type, :Instance_Type, :Instance_ID], ["AUTO", "V1", "V2", "V3"], "../plots/1_Thread_Small_Instances_Performance_Profiles.eps")

single_thread_results_of_small_instances = single_thread_results_of_small_instances[[single_thread_results_of_small_instances[:Algorithm][i] in ["AUTO", "V3", "V4", "V5", "V6", "V7", "V8", "V12", "V15", "V16"] for i in 1:nrow(single_thread_results_of_small_instances)],:]

f = plt[:figure](figsize=(20, 10))
ax = sns.boxplot(x=single_thread_results_of_small_instances[:Algorithm], y=single_thread_results_of_small_instances[:Total_Time])
ax[:set_yscale]("log")
xlabel("Algorithms")
ylabel("Computing Time (secs)")
title("Small Instances")
f[:savefig]("../plots/1_Thread_Small_Instances_Computing_Time.eps")

f = plt[:figure](figsize=(20, 10))
ax = sns.boxplot(x=single_thread_results_of_small_instances[:Algorithm], y=single_thread_results_of_small_instances[:Ratio_to_Minimum_Time])
ax[:set_yscale]("log")
xlabel("Algorithms")
ylabel("Ratio to Minimum Time")
title("Small Instances")
f[:savefig]("../plots/1_Thread_Small_Instances_Ratio_to_Minimum_Time.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results_of_small_instances[:Algorithm], y=single_thread_results_of_small_instances[:Rank_of_the_Algorithm])
xlabel("Algorithms")
ylabel("Rank of the Algorithm")
title("Small Instances")
f[:savefig]("../plots/1_Thread_Small_Instances_Rank_of_the_Algorithm.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results_of_small_instances[:Algorithm], y=single_thread_results_of_small_instances[:IP_Solved])
xlabel("Algorithms")
ylabel("IPs Solved")
title("Small Instances")
f[:savefig]("../plots/1_Thread_Small_Instances_IP_Solved.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results_of_small_instances[:Algorithm], y=single_thread_results_of_small_instances[:Percent_of_Non_Dom_Pts_Found])
xlabel("Algorithms")
ylabel("Percentage of Nondominated Points Found")
title("Small Instances")
f[:savefig]("../plots/1_Thread_Small_Instances_Percent_of_Non_Dom_Pts_Found.eps")

#####################################################################
### Plots for Large Instances                                     ###
#####################################################################

plt_performance_profiles(single_thread_results_of_large_instances, [:Problem_Type, :Instance_Type, :Instance_ID], ["AUTO", "V1", "V2", "V3"], "../plots/1_Thread_Large_Instances_Performance_Profiles.eps")

single_thread_results_of_large_instances = single_thread_results_of_large_instances[[single_thread_results_of_large_instances[:Algorithm][i] in ["AUTO", "V5", "V6", "V7", "V8", "V9", "V10", "V13", "V15", "V16"] for i in 1:nrow(single_thread_results_of_large_instances)],:]

f = plt[:figure](figsize=(20, 10))
ax = sns.boxplot(x=single_thread_results_of_large_instances[:Algorithm], y=single_thread_results_of_large_instances[:Total_Time])
ax[:set_yscale]("log")
xlabel("Algorithms")
ylabel("Computing Time (secs)")
title("Large Instances")
f[:savefig]("../plots/1_Thread_Large_Instances_Computing_Time.eps")

f = plt[:figure](figsize=(20, 10))
ax = sns.boxplot(x=single_thread_results_of_large_instances[:Algorithm], y=single_thread_results_of_large_instances[:Ratio_to_Minimum_Time])
ax[:set_yscale]("log")
xlabel("Algorithms")
ylabel("Ratio to Minimum Time")
title("Large Instances")
f[:savefig]("../plots/1_Thread_Large_Instances_Ratio_to_Minimum_Time.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results_of_large_instances[:Algorithm], y=single_thread_results_of_large_instances[:Rank_of_the_Algorithm])
xlabel("Algorithms")
ylabel("Rank of the Algorithm")
title("Large Instances")
f[:savefig]("../plots/1_Thread_Large_Instances_Rank_of_the_Algorithm.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results_of_large_instances[:Algorithm], y=single_thread_results_of_large_instances[:IP_Solved])
xlabel("Algorithms")
ylabel("IPs Solved")
title("Large Instances")
f[:savefig]("../plots/1_Thread_Large_Instances_IP_Solved.eps")

f = plt[:figure](figsize=(20, 10))
sns.boxplot(x=single_thread_results_of_large_instances[:Algorithm], y=single_thread_results_of_large_instances[:Percent_of_Non_Dom_Pts_Found])
xlabel("Algorithms")
ylabel("Percentage of Nondominated Points Found")
title("Large Instances")
f[:savefig]("../plots/1_Thread_Large_Instances_Percent_of_Non_Dom_Pts_Found.eps")


text:
module Plotting

using Plots; pyplot()
using Reproduce
using FileIO
using Statistics
using ProgressMeter
using JLD2
using LaTeXStrings
using StatsPlots

include("./plot_utils.jl")
GPU = GeneralPlotUtils
LU = LabelUtils

 default(titlefont = (20, "times"), legendfontsize = 10, guidefont = (18, :black), tickfont = 12, foreground_color_legend = nothing)

function load_data(data_home)
    ic = ItemCollection(joinpath(data_home, "data"))
    return ic
end

function load_best(ic,data_key)
    algo_divisor_keys = ["behaviour_learner", "demon_learner", "demon_opt", "demon_update","exploration_param"]
    sweep_params = ["behaviour_eta", "demon_eta", "alpha_init", "eta","emphasis_clip_threshold"]

    algo_specs_full = GPU.split_algo(ic, algo_divisor_keys)
    all_algos_ics = [search(ic,algo_spec) for algo_spec in algo_specs_full]

    valid_algos_ind = findall(a -> length(a) != 0, all_algos_ics)

    algo_specs = algo_specs_full[valid_algos_ind]
    algo_ics = all_algos_ics[valid_algos_ind]
    best_per_algo_ics = []
    for (i,algo_ic) in enumerate(algo_ics)
        push!(best_per_algo_ics, GPU.get_best_final_perf(algo_ic,sweep_params, data_key, 0.1))
    end
    return best_per_algo_ics
end

function plot_rmse_per_demon(algo_ics, data_key; inds_of_interest = nothing, save_path = "")
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ics)
    else
        deepcopy(algo_ics)[inds_of_interest]
    end

    data_per_gvf = [GPU.get_stats(GPU.load_results(ic,data_key), per_gvf = true)[1] for ic in best_per_algo_ics]
    std_per_gvf = [GPU.get_stats(GPU.load_results(ic,data_key), per_gvf = true)[2] for ic in best_per_algo_ics]
    num_gvfs = 4
    ps = []
    gvf_labels = ["Distractor" "Constant" "Drifter" "Constant"]
    num_runs = length.(best_per_algo_ics)[1]
    for gvf_ind in 1:num_gvfs
        p = plot()
        for algo_ind in 1:length(data_per_gvf)
            smooth_gvf = GPU.smooth(data_per_gvf[algo_ind][gvf_ind,:],5)
            plot_params = LU.get_params(best_per_algo_ics[algo_ind])
            plot!(p,smooth_gvf, ribbon = std_per_gvf[algo_ind][gvf_ind,:] / sqrt(num_runs), size = (500,500); plot_params...)
            plot!(xlabel="Steps", ylabel = "RMSE",  title = string(gvf_labels[gvf_ind], ""))
        end
        push!(ps,p)
    end
    allp = plot(ps..., layout=(2,2), size = (1200,1200))
    figure_path = joinpath(save_path, "Demon_RMSE_$(string(data_key)).pdf")
    savefig(figure_path)
end

function plot_TE(data_path, data_key, save_path = "")
    ic = load_data(data_path)
    algo_ic = load_best(ic, data_key)
    plot_rmse(algo_ic, data_key, save_path = save_path)
end
function plot_TE_per_demon(data_path, data_key, save_path = "")
    ic = load_data(data_path)
    algo_ic = load_best(ic, data_key)
    plot_rmse_per_demon(algo_ic, data_key, save_path = save_path)
end

function plot_goal_visitations(data_path, data_key, save_path = "")
    ic = load_data(data_path)
    algo_ic = load_best(ic, data_key)
    _plot_goal_visitations(algo_ic, data_key, save_path = save_path)
end
function plot_rmse(algo_ics, data_key; inds_of_interest = nothing, save_path = "")
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ics)
    else
        deepcopy(algo_ics)[inds_of_interest]
    end

    data = [GPU.smooth(GPU.get_stats(GPU.load_results(ic,data_key))[1],10) for ic in best_per_algo_ics]
    std = [GPU.smooth(GPU.get_stats(GPU.load_results(ic,data_key))[2],1) for ic in best_per_algo_ics]

    ylabel = "RMSE"
    title = ""
    tick_scaling = 1000
    xlabel = string("Steps", " (per $(tick_scaling))")

    step_increment = if "logger_interval" in keys(best_per_algo_ics[1].items[1].parsed_args)
        best_per_algo_ics[1].items[1].parsed_args["logger_interval"]
    else
        1000
    end
    num_samples = length(data[1])
    num_runs = length.(best_per_algo_ics)[1]
    plot_params = [LU.get_params(best_per_algo_ics[i]) for i in 1:length(best_per_algo_ics)]

    xticks=collect(1:step_increment:num_samples*step_increment) / tick_scaling
    p = plot(ylabel=ylabel, title=title, xlabel=xlabel, formatter=:plain, figsize=(800,800))
    for i in 1:length(best_per_algo_ics)
        plot!(p, xticks, data[i], ribbon = std[i]/sqrt(num_runs), legend=:none; plot_params[i]...)
    end
    figure_path = joinpath(save_path, "RMSE_$(string(data_key)).pdf")
    savefig(figure_path)
end

function _plot_goal_visitations(algo_ics, data_key; inds_of_interest = nothing, save_path="")
    best_per_algo_ics = if inds_of_interest isa Nothing
        deepcopy(algo_ics)
    else
        deepcopy(algo_ics)[inds_of_interest]
    end
    episode_lengths = [GPU.load_results(ic,:episode_length, return_type = "array") for ic in best_per_algo_ics]
    visit_orders = [GPU.load_results(ic,:visit_order, return_type = "array") for ic in best_per_algo_ics]
    max_lengths = [GPU.get_min_length(arrs)-1 for arrs in episode_lengths]
    gvf_labels = ["Distractor" "Constant" "Drifter" "Constant"]

    ps = []
    for i in 1:length(visit_orders)
        visit_perc = GPU.goal_visits_per_episode(visit_orders[i], max_lengths[i])
        tmp = [GPU.smooth(visit_perc[i,:],20) for i in 1:4]
        title = LU.get_label(best_per_algo_ics[i])[:label]

        p = plot(tmp, labels = gvf_labels, xlabel="Episode Count", ylabel="Fraction of Goal Visits", ylim=(0.0,1.0), title = title, size = (800,1000), legend=:none)
        push!(ps,p)
    end
    plot(ps...)

    figure_path = joinpath(save_path, "Goal_Visitation_$(string(data_key)).pdf")
    savefig(figure_path)
end

end

ext: jl

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"""
$(SIGNATURES)

Calculates the code phase error in chips.
"""
function dll_disc(
    system::AbstractGNSS,
    correlator,
    correlator_sample_shifts,
    early_late_index_shift,
    code_phase_delta
)
    E = abs(get_early(correlator, correlator_sample_shifts, early_late_index_shift))
    L = abs(get_late(correlator, correlator_sample_shifts, early_late_index_shift))
    distance_between_early_and_late =
        get_early_late_sample_spacing(correlator_sample_shifts, early_late_index_shift) *
        code_phase_delta
    (E - L) / (E + L) / (2 * (2 - distance_between_early_and_late))
end

"""
$(SIGNATURES)

Calculates the carrier phase error in radians.
"""
function pll_disc(system::AbstractGNSS, correlator, correlator_sample_shifts)
    p = get_prompt(correlator, correlator_sample_shifts)
    atan(imag(p) / real(p))
end


text: println("OK")

max_stars_repo_name: ShaneLee/piston

max_issues_repo_name: ShaneLee/piston

max_forks_repo_name: ShaneLee/piston


text:
function L2Frame(arg0::CelestialBody, arg1::CelestialBody)
    return L2Frame((CelestialBody, CelestialBody), arg0, arg1)
end


text:
function scale_nodes(nodes::Array{T,1},domain::Array{T,1}) where {T<:AbstractFloat}

  @inbounds for i = 1:length(nodes)
    nodes[i] = domain[2] + (1.0+nodes[i])*(domain[1]-domain[2])/2
  end

  return nodes

end

function scale_nodes(nodes::Array{T,2},domain::Array{T,2}) where {T<:AbstractFloat}

  @inbounds for i = 1:size(nodes,1)
    @inbounds for j = 1:size(nodes,2)
      nodes[i,j] = domain[2,j] + (1.0+nodes[i,j])*(domain[1,j]-domain[2,j])/2
    end
  end

  return nodes

end

function normalize_node(node::T,domain::Array{T,1}) where {T<:AbstractFloat}

  if domain[1] == domain[2]
    norm_node = zero(T)
    return norm_node
  else
    norm_node = 2.0*(node-domain[2])/(domain[1]-domain[2])-1.0
    return norm_node
  end

end

function normalize_node(node::Array{T,1},domain::Array{T,1}) where {T<:AbstractFloat}

  norm_nodes = similar(node)
  @inbounds for i in eachindex(node)
    norm_nodes[i] = normalize_node(node[i],domain)
  end

  return norm_nodes

end


text:
using Manifolds, Manopt, Test, ManifoldsBase

@testset "Test Costs" begin
    M = Sphere(2)
    N = PowerManifold(M, NestedPowerRepresentation(), 3, 3)
    f = repeat([[1.0, 0.0, 0.0]], 3, 3)
    x = repeat([[1.0, 0.0, 0.0]], 3, 3)
    x[2, 1] = [0.0, 1.0, 0.0]

    @test costIntrICTV12(N, f, f, f, 0.0, 0.0) == 0.0
    @test costIntrICTV12(N, f, x, x, 0.0, 0.0) == 1 / 2 * distance(N, x, f)^2
    @test costIntrICTV12(N, x, x, x, 2.0, 0.5) ≈ costTV2(N, x) + costTV(N, x)
    @test costIntrICTV12(N, x, x, x, 1.0, 0.0) ≈ costTV2(N, x)
    @test costIntrICTV12(N, x, x, x, 1.0, 1.0) ≈ costTV(N, x)
    @test costL2TV2(N, f, 1.0, x) == 1 / 2 * distance(N, f, x)^2 + 1.0 * costTV2(N, x)
    #
    @test costL2TV(N, f, 1.0, f) ≈ 0.0
    @test costL2TV(N, x, 1.0, x) ≈ 3 * π / 2
    @test costL2TVTV2(N, f, 0.0, 1.0, x) ≈ 1 / 2 * distance(N, x, f)^2 + costTV2(N, x)
    @test costL2TVTV2(N, f, 1.0, 1.0, x) ≈
        1 / 2 * distance(N, x, f)^2 + costTV(N, x) + costTV2(N, x)

    @test costTV2(M, Tuple(x[1:3, 1])) ≈ π / 2
    @test costTV(N, x, 1, 2) ≈ sqrt(5 / 4) * π
    @test sum(costTV2(N, x, 1, false)) == costTV2(N, x)
    @test costTV2(N, f, 2) == 0
end


text:
module test_stdlib_test

using Test # @test @inferred

f() = rand() > .5 ? 0 : 0.

@test Base.return_types(f, ()) == [Union{Float64, Int}]
@test_throws ErrorException @inferred(f())
@test_skip 1 == 2+3

@testset "one" begin
    @test true
    @testset "two" begin
        @test true
    end
end

@test length(GenericString("a")) == 1

end # module test_stdlib_test


text:
using SimpleHeatmaps
using Test

@testset "SimpleHeatmaps.jl" begin
    # Write your own tests here.
end


text:
export AbstractGenomicFeature,
    name,
    chrom,
    bpStart,
    bpEnd,
    strand,
    length,
    overlaps,
    contains,
    iscontained,
    leftof,
    rightof

"""
AbstractGenomicFeature

type implementing methods:
- name
- chrom
- bpStart
- bpEnd
- strand
"""
abstract type AbstractGenomicFeature end
# name, chrom, bpStart, bpEnd, strand

name(x::AbstractGenomicFeature) = ""
chrom(x::AbstractGenomicFeature) = x.chrom
bpStart(x::AbstractGenomicFeature) = x.bpStart
bpEnd(x::AbstractGenomicFeature) = x.bpEnd
strand(x::AbstractGenomicFeature) = x.strand

# --- abstraction boundary ----------------------------

Base.isless(x::AbstractGenomicFeature, y::AbstractGenomicFeature) =
    (chrom(x), bpStart(x), bpEnd(x)) < (chrom(y), bpStart(y), bpEnd(y))

Base.length(x::AbstractGenomicFeature)  = bpEnd(x) - bpStart(x) + 1

overlaps(x::AbstractGenomicFeature, y::AbstractGenomicFeature)  =
    (chrom(x) == chrom(y) && bpStart(x) <= bpEnd(y) && bpEnd(x) >= bpStart(y))

contains(x::AbstractGenomicFeature, y::AbstractGenomicFeature) =
    (chrom(x) == chrom(y) && bpStart(y) >= bpStart(x) && bpEnd(y) <= bpEnd(x))

iscontained(x::AbstractGenomicFeature, y::AbstractGenomicFeature) =
    contains(y, x)

leftof(x::AbstractGenomicFeature, y::AbstractGenomicFeature) =
    (chrom(x) < chrom(y))||(chrom(x) == chrom(y) && bpEnd(x) < bpStart(y))

rightof(x::AbstractGenomicFeature, y::AbstractGenomicFeature) =
    leftof(y, x)


text:
#=
These implementations were ported from the wonderful DiffLinearAlgebra
package (https://github.com/invenia/DiffLinearAlgebra.jl).
=#

using LinearAlgebra: BlasFloat

_zeros(x) = fill!(similar(x), zero(eltype(x)))

#####
##### `BLAS.dot`
#####

frule((Δself, Δx, Δy), ::typeof(BLAS.dot), x, y) = frule((Δself, Δx, Δy), dot, x, y)

rrule(::typeof(BLAS.dot), x, y) = rrule(dot, x, y)

function rrule(::typeof(BLAS.dot), n, X, incx, Y, incy)
    Ω = BLAS.dot(n, X, incx, Y, incy)
    function blas_dot_pullback(Ω̄)
        ΔΩ = unthunk(Ω̄)
        if ΔΩ isa ZeroTangent
            ∂X = ZeroTangent()
            ∂Y = ZeroTangent()
        else
            ∂X = @thunk scal!(n, ΔΩ, blascopy!(n, Y, incy, _zeros(X), incx), incx)
            ∂Y = @thunk scal!(n, ΔΩ, blascopy!(n, X, incx, _zeros(Y), incy), incy)
        end
        return (NoTangent(), NoTangent(), ∂X, NoTangent(), ∂Y, NoTangent())
    end
    return Ω, blas_dot_pullback
end

#####
##### `BLAS.nrm2`
#####

function frule((_, Δx), ::typeof(BLAS.nrm2), x)
    Ω = BLAS.nrm2(x)
    s = ifelse(iszero(Ω), one(Ω), Ω)
    ∂Ω = if x isa Real
        BLAS.dot(x, Δx) / s
    else
        sum(y -> _realconjtimes(y...), zip(x, Δx)) / s
    end
    return Ω, ∂Ω
end

function rrule(::typeof(BLAS.nrm2), x)
    Ω = BLAS.nrm2(x)
    function nrm2_pullback(ΔΩ)
        return NoTangent(), x .* (real(ΔΩ) / ifelse(iszero(Ω), one(Ω), Ω))
    end
    return Ω, nrm2_pullback
end

function rrule(::typeof(BLAS.nrm2), n, X, incx)
    Ω = BLAS.nrm2(n, X, incx)
    nrm2_pullback(::ZeroTangent) = (NoTangent(), NoTangent(), ZeroTangent(), NoTangent())
    function nrm2_pullback(ΔΩ)
        # BLAS.scal! requires s has the same eltype as X
        s = eltype(X)(real(ΔΩ) / ifelse(iszero(Ω), one(Ω), Ω))
        ∂X = scal!(n, s, blascopy!(n, X, incx, _zeros(X), incx), incx)
        return (NoTangent(), NoTangent(), ∂X, NoTangent())
    end
    return Ω, nrm2_pullback
end

#####
##### `BLAS.asum`
#####

function frule((_, Δx), ::typeof(BLAS.asum), x)
    ∂Ω = sum(zip(x, Δx)) do (xi, Δxi)
        return _realconjtimes(_signcomp(xi), Δxi)
    end
    return BLAS.asum(x), ∂Ω
end

function rrule(::typeof(BLAS.asum), x)
    function asum_pullback(ΔΩ)
        return (NoTangent(), _signcomp.(x) .* real(ΔΩ))
    end
    return BLAS.asum(x), asum_pullback
end

function rrule(::typeof(BLAS.asum), n, X, incx)
    Ω = BLAS.asum(n, X, incx)
    asum_pullback(::ZeroTangent) = (NoTangent(), NoTangent(), ZeroTangent(), NoTangent())
    function asum_pullback(ΔΩ)
        # BLAS.scal! requires s has the same eltype as X
        s = eltype(X)(real(ΔΩ))
        ∂X = scal!(n, s, blascopy!(n, _signcomp.(X), incx, _zeros(X), incx), incx)
        return (NoTangent(), NoTangent(), ∂X, NoTangent())
    end
    return Ω, asum_pullback
end

# component-wise sign, e.g. sign(x) + i sign(y)
@inline _signcomp(x::Real) = sign(x)
@inline _signcomp(x::Complex) = Complex(sign(real(x)), sign(imag(x)))

#####
##### `BLAS.gemv`
#####

function rrule(::typeof(gemv), tA::Char, α::T, A::AbstractMatrix{T},
               x::AbstractVector{T}) where T<:BlasFloat
    y = gemv(tA, α, A, x)
    function gemv_pullback(Ȳ)
        ȳ = unthunk(Ȳ)
        if uppercase(tA) === 'N'
            ∂A = InplaceableThunk(
                @thunk(α' * ȳ * x'),
                Ā -> ger!(α', ȳ, x, Ā)
            )
            ∂x = InplaceableThunk(
                @thunk(gemv('C', α', A, ȳ)),
                x̄ -> gemv!('C', α', A, ȳ, one(T), x̄)
            )
        elseif uppercase(tA) === 'C'
            ∂A = InplaceableThunk(
                @thunk(α * x * ȳ'),
                Ā -> ger!(α, x, ȳ, Ā)
            )
            ∂x = InplaceableThunk(
                @thunk(gemv('N', α', A, ȳ)),
                x̄ -> gemv!('N', α', A, ȳ, one(T), x̄)
            )
        else  # uppercase(tA) === 'T'
            ∂A = InplaceableThunk(
                @thunk(conj(α * x * ȳ')),
                Ā -> conj!(ger!(α, x, ȳ, conj!(Ā)))
            )
            ∂x = InplaceableThunk(
                @thunk(gemv('N', α', conj(A), ȳ)),
                x̄ -> gemv!('N', α', conj(A), ȳ, one(T), x̄)
            )
        end
        return (NoTangent(), NoTangent(), @thunk(dot(y, ȳ) / α'), ∂A, ∂x)
    end
    return y, gemv_pullback
end

function rrule(
    ::typeof(gemv), tA::Char, A::AbstractMatrix{T}, x::AbstractVector{T}
) where T<:BlasFloat
    y, inner_pullback = rrule(gemv, tA, one(T), A, x)
    function gemv_pullback(Ȳ)
        (_, dtA, _, dA, dx) = inner_pullback(Ȳ)
        return (NoTangent(), dtA, dA, dx)
    end
    return y, gemv_pullback
end

#####
##### `BLAS.gemm`
#####

function rrule(
    ::typeof(gemm), tA::Char, tB::Char, α::T, A::AbstractMatrix{T}, B::AbstractMatrix{T}
) where T<:BlasFloat
    C = gemm(tA, tB, α, A, B)
    function gemm_pullback(Cbar)
        C̄ = unthunk(Cbar)
        β = one(T)
        if uppercase(tA) === 'N'
            if uppercase(tB) === 'N'
                ∂A = InplaceableThunk(
                    @thunk(gemm('N', 'C', α', C̄, B)),
                    Ā -> gemm!('N', 'C', α', C̄, B, β, Ā)
                )
                ∂B = InplaceableThunk(
                    @thunk(gemm('C', 'N', α', A, C̄)),
                    B̄ -> gemm!('C', 'N', α', A, C̄, β, B̄)
                )
            elseif uppercase(tB) === 'C'
                ∂A = InplaceableThunk(
                    @thunk(gemm('N', 'N', α', C̄, B)),
                    Ā -> gemm!('N', 'N', α', C̄, B, β, Ā)
                )
                ∂B = InplaceableThunk(
                    @thunk(gemm('C', 'N', α, C̄, A)),
                    B̄ -> gemm!('C', 'N', α, C̄, A, β, B̄)
                )
            else  # uppercase(tB) === 'T'
                ∂A = InplaceableThunk(
                    @thunk(gemm('N', 'N', α', C̄, conj(B))),
                    Ā -> gemm!('N', 'N', α', C̄, conj(B), β, Ā)
                )
                ∂B = InplaceableThunk(
                    @thunk(conj(gemm('C', 'N', α, C̄, A))),
                    B̄ -> conj!(gemm!('C', 'N', α, C̄, A, β, conj!(B̄)))
                )
            end
        elseif uppercase(tA) === 'C'
            if uppercase(tB) === 'N'
                ∂A = InplaceableThunk(
                    @thunk(gemm('N', 'C', α, B, C̄)),
                    Ā -> gemm!('N', 'C', α, B, C̄, β, Ā)
                )
                ∂B = InplaceableThunk(
                    @thunk(gemm('N', 'N', α', A, C̄)),
                    B̄ -> gemm!('N', 'N', α', A, C̄, β, B̄)
                )
            elseif uppercase(tB) === 'C'
                ∂A = InplaceableThunk(
                    @thunk(gemm('C', 'C', α, B, C̄)),
                    Ā -> gemm!('C', 'C', α, B, C̄, β, Ā)
                )
                ∂B = InplaceableThunk(
                    @thunk(gemm('C', 'C', α, C̄, A)),
                    B̄ -> gemm!('C', 'C', α, C̄, A, β, B̄)
                )
            else  # uppercase(tB) === 'T'
                ∂A = InplaceableThunk(
                    @thunk(gemm('T', 'C', α, B, C̄)),
                    Ā -> gemm!('T', 'C', α, B, C̄, β, Ā)
                )
                ∂B = InplaceableThunk(
                    @thunk(gemm('T', 'T', α', C̄, A)),
                    B̄ -> gemm!('T', 'T', α', C̄, A, β, B̄)
                )
            end
        else  # uppercase(tA) === 'T'
            if uppercase(tB) === 'N'
                ∂A = InplaceableThunk(
                    @thunk(conj(gemm('N', 'C', α, B, C̄))),
                    Ā -> conj!(gemm!('N', 'C', α, B, C̄, β, conj!(Ā)))
                )
                ∂B = InplaceableThunk(
                    @thunk(gemm('N', 'N', α', conj(A), C̄)),
                    B̄ -> gemm!('N', 'N', α', conj(A), C̄, β, B̄)
                )
            elseif uppercase(tB) === 'C'
                ∂A = InplaceableThunk(
                    @thunk(gemm('T', 'T', α', B, C̄)),
                    Ā -> gemm!('T', 'T', α', B, C̄, β, Ā)
                )
                ∂B = InplaceableThunk(
                    @thunk(gemm('C', 'T', α, C̄, A)),
                    B̄ -> gemm!('C', 'T', α, C̄, A, β, B̄)
                )
            else  # uppercase(tB) === 'T'
                ∂A = InplaceableThunk(
                    @thunk(gemm('C', 'T', α', B, C̄)),
                    Ā -> gemm!('C', 'T', α', B, C̄, β, Ā)
                )
                ∂B = InplaceableThunk(
                    @thunk(gemm('T', 'C', α', C̄, A)),
                    B̄ -> gemm!('T', 'C', α', C̄, A, β, B̄)
                )
            end
        end
        return (NoTangent(), NoTangent(), NoTangent(), @thunk(dot(C, C̄) / α'), ∂A, ∂B)
    end
    return C, gemm_pullback
end

function rrule(
    ::typeof(gemm), tA::Char, tB::Char, A::AbstractMatrix{T}, B::AbstractMatrix{T}
) where T<:BlasFloat
    C, inner_pullback = rrule(gemm, tA, tB, one(T), A, B)
    function gemm_pullback(Ȳ)
        (_, dtA, dtB, _, dA, dB) = inner_pullback(Ȳ)
        return (NoTangent(), dtA, dtB, dA, dB)
    end
    return C, gemm_pullback
end


text:
using Test
using FindDefinition

include("bars.jl")

barspath = joinpath(@__DIR__, "bars.jl")

@testset "single macro call" begin
    bar1_methoddefs = FindDefinition.find_definitions(only(methods(Bar1.bar)))
    bar1_func_lnns = finddefs(Bar1.bar)
    @test length(bar1_methoddefs) == 1
    @test [d.lnn for d in bar1_methoddefs] == bar1_func_lnns
    lnn = only(bar1_func_lnns)
    @test string(lnn.file) == barspath
    @test lnn.line == 6
end

@testset "several calls" begin
    meths = methods(Bar2.bar)
    bar2_methoddefs = FindDefinition.find_definitions.(meths)
    bar2_func_lnns = finddefs(Bar2.bar)
    @test length(bar2_methoddefs) == 2
    @test length(bar2_methoddefs[1]) == 2
    @test length(bar2_methoddefs[2]) == 1
    @test length(bar2_func_lnns) == 2
    @test all( all( string(d.lnn.file) == barspath for d in defs )
                for defs in bar2_methoddefs )
    @test bar2_methoddefs[1][1].lnn.line == 12
    @test bar2_methoddefs[1][2].lnn.line == 13
    @test bar2_methoddefs[2][1].lnn.line == 14

    @test bar2_func_lnns[1].line == 13
    @test bar2_func_lnns[2].line == 14
end

@testset "several definitions in one macro" begin
    meths = methods(Bar2.multibar)
    @assert length(meths) == 3
    defs = finddefs(Bar2.multibar)
    @test length(defs) == 3
    @test defs[1] == defs[2] == defs[3]
    @test string(defs[1].file) == barspath
    @test defs[1].line == 15
end


text:
# ---
# title: 529. Minesweeper
# id: problem529
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/minesweeper/description/>
# hidden: true
# ---
#
# Let's play the minesweeper game
# ([Wikipedia](https://en.wikipedia.org/wiki/Minesweeper_\(video_game\)),
# [online game](http://minesweeperonline.com))!
#
# You are given a 2D char matrix representing the game board. **' M'**
# represents an **unrevealed** mine, **' E'** represents an **unrevealed** empty
# square, **' B'** represents a **revealed** blank square that has no adjacent
# (above, below, left, right, and all 4 diagonals) mines, **digit** ('1' to '8')
# represents how many mines are adjacent to this **revealed** square, and
# finally **' X'** represents a **revealed** mine.
#
# Now given the next click position (row and column indices) among all the
# **unrevealed** squares ('M' or 'E'), return the board after revealing this
# position according to the following rules:
#
#   1. If a mine ('M') is revealed, then the game is over - change it to **' X'**.
#   2. If an empty square ('E') with **no adjacent mines** is revealed, then change it to revealed blank ('B') and all of its adjacent **unrevealed** squares should be revealed recursively.
#   3. If an empty square ('E') with **at least one adjacent mine** is revealed, then change it to a digit ('1' to '8') representing the number of adjacent mines.
#   4. Return the board when no more squares will be revealed.
#
#
#
# **Example 1:**
#
#
#
#     Input:
#
#     [['E', 'E', 'E', 'E', 'E'],
#      ['E', 'E', 'M', 'E', 'E'],
#      ['E', 'E', 'E', 'E', 'E'],
#      ['E', 'E', 'E', 'E', 'E']]
#
#     Click : [3,0]
#
#     Output:
#
#     [['B', '1', 'E', '1', 'B'],
#      ['B', '1', 'M', '1', 'B'],
#      ['B', '1', '1', '1', 'B'],
#      ['B', 'B', 'B', 'B', 'B']]
#
#     Explanation:
#     ![](https://assets.leetcode.com/uploads/2018/10/12/minesweeper_example_1.png)
#
#
# **Example 2:**
#
#
#
#     Input:
#
#     [['B', '1', 'E', '1', 'B'],
#      ['B', '1', 'M', '1', 'B'],
#      ['B', '1', '1', '1', 'B'],
#      ['B', 'B', 'B', 'B', 'B']]
#
#     Click : [1,2]
#
#     Output:
#
#     [['B', '1', 'E', '1', 'B'],
#      ['B', '1', 'X', '1', 'B'],
#      ['B', '1', '1', '1', 'B'],
#      ['B', 'B', 'B', 'B', 'B']]
#
#     Explanation:
#     ![](https://assets.leetcode.com/uploads/2018/10/12/minesweeper_example_2.png)
#
#
#
#
# **Note:**
#
#   1. The range of the input matrix's height and width is [1,50].
#   2. The click position will only be an unrevealed square ('M' or 'E'), which also means the input board contains at least one clickable square.
#   3. The input board won't be a stage when game is over (some mines have been revealed).
#   4. For simplicity, not mentioned rules should be ignored in this problem. For example, you **don 't** need to reveal all the unrevealed mines when the game is over, consider any cases that you will win the game or flag any squares.
#
#
## @lc code=start
using LeetCode

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


text:
using YAActL, .Iterators

@enum Q s t u v w   # describe states

# implement the transition function/methods
δ(::Val{s}, ::Val{'a'}) = (t,v,w)
δ(::Val{s}, ::Any)      = (s,)
δ(::Val{t}, ::Val{'a'}) = (u,v)
δ(::Val{t}, ::Val{'b'}) = (u,v)
δ(::Val{t}, ::Any)      = (t,)
δ(::Val{u}, ::Val{'a'}) = (s,w)
δ(::Val{u}, ::Val{'b'}) = (s,w)
δ(::Val{u}, ::Any)      = (u,)
δ(::Val{v}, ::Val{'b'}) = (t,v)
δ(::Val{v}, ::Val{'c'}) = (t,v)
δ(::Val{v}, ::Any)      = (v,)
δ(::Val{w}, ::Val{'b'}) = (t,v)
δ(::Val{w}, ::Any)      = (w,)
δ(q::Q, c::Char) = δ(Val(q), Val(c))
δ(qs::Tuple, c::Char) = [δ(q,c) for q in qs]|>flatten|>Set|>Tuple

# this is a simple iteration over state
function check(str::String)
    qs = s
    for c in str
        qs = δ(qs, c)
    end
    intersect(qs,(u,v)) |> !isempty
end

# this works with the actor in state dispatch
function check(lk::Link, str::String)
    update!(mynfa, s)
    foreach(c->cast!(lk,c), str)
    intersect(query!(lk), (u,v)) |> !isempty
end

mynfa = Actor(δ)
set!(mynfa, state)

check(mynfa, "aabc")

check(mynfa, "bbb")


text:
using eFEMpart

# test fem_generate_map
rect = [-1.0,1.0,-1.0,1.0]
N    = 64
mesh = squareMesh(rect,N,2)
L    = 0.1
totalBounds = [rect[1]-L,rect[2]+L,rect[3]-L,rect[4]+L]
femCLmap = StokesParticles.femGenerateMap(mesh,totalBounds,L)

# test getElementFromNode!
rect = [-1.0,1.0,-1.0,1.0]
N    = 5
mesh = squareMesh(rect,N,1)
elementArray = zeros(Int,4)
nodeInd = 12
@test StokesParticles.getElementFromNode!(elementArray,mesh,nodeInd) != π     # just to make sure function runs

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
using Documenter, HyperDualMatrixTools, LinearAlgebra

makedocs(
    sitename="HyperDualMatrixTools Documentation",
    # options
    modules = [HyperDualMatrixTools]
)

deploydocs(
    repo = "github.com/briochemc/HyperDualMatrixTools.jl.git",
)


text:
using DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
using Test, Distributions, SteadyStateDiffEq
println("One parameter case")
f1 = @ode_def begin
    dx = a*x - x*y
    dy = -3y + x*y
end a
u0 = [1.0,1.0]
tspan = (0.0,10.0)
prob1 = ODEProblem(f1,u0,tspan,[1.5])
sol = solve(prob1,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [Normal(1.5,0.01)]

bayesian_result = turing_inference(prob1,Tsit5(),t,data,priors;num_samples=500,
                                   syms=[:a])

@show bayesian_result

@test mean(get(bayesian_result,:a)[1]) ≈ 1.5 atol=3e-1

bayesian_result = turing_inference(prob1,Rosenbrock23(autodiff=false),t,data,priors;num_samples=500,
                                   syms=[:a])

bayesian_result = turing_inference(prob1,Rosenbrock23(),t,data,priors;num_samples=500,
                                   syms=[:a])

priors = [Normal(1.,0.01),Normal(1.,0.01),Normal(1.5,0.01)]
bayesian_result = turing_inference(prob1,Tsit5(),t,data,priors;num_samples=500,sample_u0 =true,
                                   syms=[:u1,:u2,:a])

@test mean(get(bayesian_result,:a)[1]) ≈ 1.5 atol=3e-1
@test mean(get(bayesian_result,:u1)[1]) ≈ 1.0 atol=3e-1
@test mean(get(bayesian_result,:u2)[1]) ≈ 1.0 atol=3e-1

sol = solve(prob1,Tsit5(),save_idxs=[1])
randomized = VectorOfArray([(sol(t[i]) + .01 * randn(1)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [Normal(1.5,0.01)]
bayesian_result = turing_inference(prob1,Tsit5(),t,data,priors;num_samples=500,
                                   syms=[:a],save_idxs=[1])

@test mean(get(bayesian_result,:a)[1]) ≈ 1.5 atol=3e-1

priors = [Normal(1.,0.01),Normal(1.5,0.01)]
bayesian_result = turing_inference(prob1,Tsit5(),t,data,priors;num_samples=500,sample_u0 =true,
                                   syms=[:u1,:a],save_idxs=[1])

@test mean(get(bayesian_result,:a)[1]) ≈ 1.5 atol=3e-1
@test mean(get(bayesian_result,:u1)[1]) ≈ 1.0 atol=3e-1

println("Four parameter case")
f2 = @ode_def begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob2 = ODEProblem(f2,u0,tspan,p)
sol = solve(prob2,Tsit5())
t = collect(range(1,stop=10,length=10))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
priors = [truncated(Normal(1.5,0.01),0,2),truncated(Normal(1.0,0.01),0,1.5),
          truncated(Normal(3.0,0.01),0,4),truncated(Normal(1.0,0.01),0,2)]

bayesian_result = turing_inference(prob2,Tsit5(),t,data,priors;num_samples=500,
                                   syms = [:a,:b,:c,:d])

@show bayesian_result

@test mean(get(bayesian_result,:a)[1]) ≈ 1.5 atol=3e-1
@test mean(get(bayesian_result,:b)[1]) ≈ 1.0 atol=3e-1
@test mean(get(bayesian_result,:c)[1]) ≈ 3.0 atol=3e-1
@test mean(get(bayesian_result,:d)[1]) ≈ 1.0 atol=3e-1

println("Steady state problem")
function f(du,u,p,t)
    α = p[1]
    du[1] = 2 -  α*u[1]
    du[2] = u[1] - 4u[2]
end

p = [2.0]
u0 = zeros(2)
s_prob = SteadyStateProblem(f,u0,p)
s_sol = solve(s_prob,SSRootfind())
s_sol = solve(s_prob,DynamicSS(Tsit5(),abstol=1e-4,reltol=1e-3))

# true data is 1.00, 0.25
data = [1.05, 0.23]
priors = [truncated(Normal(2.0,0.2),0,3)]
bayesian_result = turing_inference(s_prob,DynamicSS(Tsit5(),abstol=1e-4,reltol=1e-3),
                                   nothing,data,priors;
                                   num_samples=500,
                                   maxiters = 1e6,
                                   syms = [:α])

@test mean(get(bayesian_result,:α)[1]) ≈ 2.0 atol=3e-1


text:
# This script is a (currently unsuccessful) attempt to implement a partial differential
# equation describing changes in concentration with respect to longitudinal distance along
# a stream channel. Changes in concentration result from both advection and a first-order
# reaction in the hyporheic zone, which exchanges with the channel in such a way that
# the ages of water returning to the channel are distributed as a power law.

##### ----- Define Parameters, Initial Conditions, Boundary Conditions ----- #####
τ₋ = 1.0 # minimum residence time in seconds
τ₊ = 1000.0 # maximum residence time in seconds
α = 1.6 # power-law exponent
k = 0.0 # first-order rate constant in seconds^-1
S = 1.0 # volume of water in hyporheic zone
# Calculate rate of hyporheic exchange (volume per time)
q = qCalc_powerLaw(S, τ₋, τ₊, α)
v = 1 # Velocity
p = [τ₋, τ₊, α, k, q, v] # parameter vector
tspan = (0.0, 10.0) # simulation timespan
xₙ = 10 # physical length of simulation
N = 999 # number of interior points
Δx = xₙ/(N+1) # grid spacing along x dimension (longitudinal)
x = range(0, step = Δx, length = (N+2)) # x-points of solution, including boundaries
h(p,t) = zeros(N) # history function describing solute concentrations before the start of the simulation
C₀ = zeros(N) # initial condition
for i in 1:N
    C₀[i] = 1.0 - i*Δx/xₙ
end

##### ----- Finite Difference Operator Method ----- #####
ΔUpwind = DiffEqOperators.UpwindDifference(1,4,Δx,xₙ) # first derivative operator, fourth-order approximation
Dirichletbc = DiffEqOperators.DirichletBC(1.0,0.0) # boundary condition operator

@advectReact(ΔUpwind,Dirichletbc)

# Build the DDE problem
prob = DiffEqBase.ODEProblem(advectionReaction!, C₀, tspan, p)

# Solve problem
solution = solve(prob,Rosenbrock23())

# Plot results
t = solution.t
u = solution.u

plt = plot()
for i in 1:length(t)
    plot!(x, Dirichletbc*u[[i]][1] ,label="t=$(t[i])")
end
display(plt)

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
##
## --------------------------------------------------
##
## Classic Continuos Orthogonal Polynomials
abstract type AbstractCCOP{T,N} <: AbstractCOP{T,N} end
"""
    AbstractCCOP{T,N} <:  AbstractCOP{T,N}

Following [Koepf and Schmersau](https://arxiv.org/pdf/math/9703217.pdf), a family `y(x)=p_n(x)=k_x⋅x^n +  ...`
for  `n  ∈  {0, 1,…}, k_n ≠ 0` of polynomials is a family of classic *continuous* orthogonal polynomials if each is  a
solution of the differential equation

`(a⋅x²+b⋅x+c) ⋅ yᵢ'' + (d⋅x + e) ⋅ yᵢ' + λᵢ⋅ yᵢ = 0.`

A family is characterized, up to choice of leading term, by the 5 coefficients: `a,b,c,d,e`.
Let `σ = (a⋅x²+b⋅x+c)`, `τ = (d⋅x + e)`.

From these  5  coefficients several structural  equations are represented. For example
the three-point recusion.

`P₍ᵢ₊₁) = (Aᵢ⋅x + Bᵢ) * Pᵢ - Cᵢ *  P₍ᵢ₋₁₎`,

where `Aᵢ,Bᵢ,Cᵢ` can be represented in formulas involving just  `a,b,c,d,e` and `i`.

Rearranging   gives the structural equation:

`x⋅p_n   = [an, bn, cn] ⋅ [p_{n+1}, p_n, p_{n-1}]`  (Eqn (7))

The other structural equations are (equation  references are from Koepf and Schmersau):

`σ⋅p'_n  = [αn, βn, γn] ⋅  [p_{n+1}, p_n, p_{n-1}]` (Eqn (9), n ≥ 1)

`p_n = [ân, b̂n, ĉn]  ⋅  [p'_{n+1}, p'_n, p'_{n-1}]` (Eqn (19))

`x⋅p'_n  = [αᴵn, βᴵn, γᴵn] ⋅  [p'_{n+1}, p'_n, p'_{n-1}]` (Eqn  (14))

Using (7), Clenshaw polynomial evaluation using the three  point recursion is defined.

Using (19), expressions for derivatives are found.

Using  (19), expressions for integration are found (p7).

Using their theorems 2,4, and 5, connection coefficients, `C(n,m)` satisfying
`P_n(x) =  ∑  C(n,m)  Q_m(x) (n ≥ 0, 0 ≤  m ≤ n)` are  found. These
allow  fallback  definitions for `convert(Polynomial,p)`,  `convert(P, p::Polynomial)`,
`convert(P{α…}, p::P(β…))` and through composition polynomial  multiplication,  `p*q`.

Subtypes of `AbstractCCOP` are  created through  the `@register0` or  `@registerN` macros, where the  `N`  macro  is used  if parameters are  needed to describe the family.

If non-monic versions are desired, then the  leading  term can be  specified through `kn()` (which by default is defined by the  method `k1k0(P,i)`, the ratio of  `kᵢ₊₁/kᵢ`).  The `@register_monic` macro is useful  for creating  monic versions through  method delegation from the common non-monic systems. Similarly, the `@register_shifted` macro is useful  to provide shifted versions (cf. [`ShiftedLegendre`](@ref)).

Registering a system, defining an `abcde` method, and optionally
defining `k1k0` is usually sufficient to define a new system, though
the general equations may need specializations when algebraic
cancellation is required.

The defaults for evaluation and multplication are generally an order
of magnitude slower than a directly defined function. For some
families this is done (e.g. `Chebyshev`,`ChebyshevU`, `Hermite`, `Laguerre`), but not all.

## Example

For this example, the value of `Bn` at `0` needs help:

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> const SP=SpecialPolynomials
SpecialPolynomials

julia> SP.@register0 MonicLegendre′ SP.AbstractCCOP0

julia> SP.:ϟ(::Type{<:MonicLegendre′}) = Legendre

julia> SP.@register_monic MonicLegendre′  # use  ϟ to delegate methods

julia> 𝐐  =  Rational{Int}
Rational{Int64}

julia> x = variable(Polynomial{𝐐})
Polynomial(x)

julia> [basis(MonicLegendre′{𝐐}, i)(x) for i  in 0:5]
6-element Array{Polynomial{Rational{Int64}},1}:
 Polynomial(1//1)
 Polynomial(x)
 Polynomial(-1//3 + x^2)
 Polynomial(-3//5*x + x^3)
 Polynomial(3//35 - 6//7*x^2 + x^4)
 Polynomial(5//21*x - 10//9*x^3 + x^5)
```

[Koekoek and Swarttouw](https://arxiv.org/pdf/math/9602214.pdf)
present an encyclopedia of formula characterizing families of
orthogonal polynomials.

"""
AbstractCCOP

# subtypes  to keep track of number of parameters
# passed to  @registerN macros
abstract type AbstractCCOP0{T,N} <: AbstractCCOP{T,N} end
abstract type AbstractCCOP1{α,T,N} <: AbstractCCOP{T,N} end
abstract type AbstractCCOP2{α,β,T,N} <: AbstractCCOP{T,N}  end
abstract type AbstractCCOP3{α,β,γ,T,N} <: AbstractCCOP{T,N}  end

# We want to  be able to strip  off T,N  or α,...,T,N
# * constructorof(P{α,..,T,N}) =  P
# * ⟒(P(α,...,T,N)) =  P(α...)  #  \upin[tab]
⟒(P::Type{<:AbstractCCOP1{α}}) where {α} = constructorof(P){α}
⟒(P::Type{<:AbstractCCOP2{α,β}}) where {α, β} = constructorof(P){α, β}
⟒(P::Type{<:AbstractCCOP3{α,β,γ}}) where {α, β,γ} = constructorof(P){α, β, γ}

# for conversion to base case
# \upstigma[tab]
ϛ(P::Type{<:AbstractCCOP}) = Polynomial

## Display
# show parameters in constructor's name
function Base.show(io::IO, mimetype::MIME"text/plain", p::P) where {α,P<:AbstractCCOP1{α}}
    print(io,"$(P.name){$(α)}(")
    printpoly(io, p, mimetype)
    print(io,")")
end
function Base.show(io::IO, mimetype::MIME"text/plain", p::P) where {α,β,P<:AbstractCCOP2{α,β}}
    print(io,"$(P.name){$(α),$(β)}(")
    printpoly(io, p, mimetype)
    print(io,")")
end

##
## -----
##
##
## -----
##

## abcde for each family
## Families       a   b   c   d       e
## Hermite        1   0   0  -2       0
## Hermiteₑ       0   0   1  -1       0
## ChebyshevT    -1   0   1  -1       0
## ChebyshevU    -1   0   1  -3       0
## Laguerre{α}    0   1   0  -1      α+1
## Bessel{α}      1   0   0  α        2
## Gegenbaeur{α} -1   0   1 -2(α+1)   0
## Jacobi{α,β}   -1   0   1 -(α+β+2) β-α
## standard

##
## Structural  equations
##
## Could move the AbstractCOP terms into `cop.jl`...

# An, Bn,  Cn
# p_{n+1} = (An*x + Bn)⋅p_n + Cn⋅p_{n-1}
An(P::Type{<:AbstractCOP}, n::Int) = Ãn(P,n) * k1k0(P, n)
function  Ãn(P::Type{<:AbstractCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    Ãn(P, a,b,c,d,e ,n)
end

Bn(P::Type{<:AbstractCOP}, n::Int) = B̃n(P,n) * k1k0(P,n)
function B̃n(P::Type{<:AbstractCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    B̃n(P, a,b,c,d,e ,n)
end

Cn(P::Type{<:AbstractCOP}, n::Int) = C̃n(P,n) * k1k_1(P, n)
function C̃n(P::Type{<:AbstractCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    C̃n(P, a,b,c,d,e,n)
end

function Ãn(P::Type{<:AbstractCCOP}, a,b,c,d,e, n::Int)
    one(eltype(P))
end

function B̃n(P::Type{<:AbstractCCOP}, a,b,c,d,e, n::Int)

    S = eltype(P)

    num = (2b*n*(a*n+d-a)-e*(-d+2a))
    den = (d+2a*n) * (d-2a+2a*n)

    iszero(den) && return B̃n(P, Val(n))

    val = (one(S) * num) / den

    val
end

function C̃n(P::Type{<:AbstractCCOP}, a,b,c,d,e, n::Int)

    S = eltype(P)

    numa = (a*n+d-2a) * n * (4c*a-b^2) + 4a^2*c -a*b^2 + a*e^2 - 4*a*c*d
    numa += d*b^2 - b*e*d + d^2*c
    num = -numa * (a*n + d - 2a) * n
    den = (d - 2a + 2a*n)^2 * (2a*n - 3a + d) * (2a*n - a + d)

    iszero(den) && return C̃n(P, Val(n))

    val = (one(S) * num)  / den

    val
end

# an, bn, cn
# x⋅pn = [an,bn,cn] ⋅ [p_{n+1},p_n,p_{n-1}]
function an(P::Type{<:AbstractCOP}, n::Int)
    1/An(P,n)
end

function bn(P::Type{<:AbstractCOP}, n::Int)
    -Bn(P,n)/An(P,n)
end

function cn(P::Type{<:AbstractCOP}, n::Int)
    Cn(P,n)/An(P,n)
end

# αn, βn,γn
# σ⋅pn' = [αn, βn,γn] ⋅ [p_{n+1},p_n,p_{n-1}]
αn(P::Type{<:AbstractCOP}, n::Int) = α̃n(P,n) / k1k0(P,n)
function α̃n(P::Type{<:AbstractCCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    val = (one(S) *  a) * n

    return val
end

βn(P::Type{<:AbstractCOP}, n::Int) = β̃n(P, n)
function β̃n(P::Type{<:AbstractCCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = -n*(a*n+d-a)*(2e*a-d*b)
    den = (d+2*a*n)*(d-2a+2a*n)

    iszero(den) &&  return β̃n(P, Val(n))

    val = (one(S) * num)  / den

    return val
end

γn(P::Type{<:AbstractCOP}, n::Int) = γ̃n(P,n) *  k1k0(P,n-1)
function γ̃n(P::Type{<:AbstractCCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = ((n-1) * (a*n + d - a) * (4c*a - b^2) + a*e^2 + d^2*c - b*e*d)
    num *= (a*n + d - a) * (a*n + d - 2a) * n
    den = (d - 2a + 2a*n)^2 * (2a*n - 3a + d) * (2a*n - a + d)
    iszero(den) &&  return γ̃n(P, Val(n))

    val = (one(S) * num) / den

    return val
end

# for integration formulas
# ân, b̂n, ĉn
# pn = [ân, b̂n, ĉn] ⋅ [p'_{n+1},p'_n,p'_{n-1}]
ân(P::Type{<:AbstractCOP}, n::Int) = ẫn(P, n) / k1k0(P,n)
function ẫn(P::Type{<:AbstractCCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    val = one(S)
    val /= n+1
    return val
end

b̂n(P::Type{<:AbstractCOP}, n::Int) =  b̂̃n(P,n)
function b̂̃n(P::Type{<:AbstractCCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = (2e*a - d*b)
    den = (d+2a*n)*(d-2a+2a*n)
    iszero(den) && return b̂̃n(P, Val(n))

    val = one(S) * num/ den
    return val

end

ĉn(P::Type{<:AbstractCOP}, n::Int) = ĉ̃n(P,n) * k1k0(P,n-1)
function ĉ̃n(P::Type{<:AbstractCCOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = ((n-1) * (a*n  + d - a) * (4c*a - b^2) + a*e^2 +  d^2*c - b*e*d) * a* n
    den = (d  - 2a + 2a*n) * (d - 2a + 2a*n) * (2a*n - 3a + d) * (2a*n - a  + d)
    iszero(den)  &&  return ĉ̃n(P, Val(n))

    val = (one(S)  *  num) / den
    return val
end

"""
    abcdeᴵ

If P_n is a CCOP, then

x pᵢ = αᵢ p₍ᵢ₊₁₎ + βᵢ pᵢ + γᵢ p₍ᵢ₋₁₎.

Sp will pᵢ' with this  choice of a,b,d,e derived from  those of  pᵢ. (Eqn  (13)
"""
function abcdeᴵ(P::Type{<:AbstractCCOP})
    a,b,c,d,e = abcde(P).a, abcde(P).b, abcde(P).c, abcde(P).d, abcde(P).e
    NamedTuple{(:a,:b,:c,:d,:e)}((a,b,c,d+2a,e+b))
end

# αᴵn, βᴵn, γᴵn  (α^*,...)
# x⋅pn' = [αᴵn, βᴵn, γᴵn] ⋅  [p'_{n+1},p'_n,p'_{n-1}]
αᴵn(P::Type{<:AbstractCOP}, n::Int) = α̃ᴵn(P,n) / k1k0(P,n)
function α̃ᴵn(P::Type{<:AbstractCOP}, n::Int)
    n  = n - 1
    a,b,c,d,e = abcdeᴵ(P)
    Aᴵn = Ãn(P,a,b,c,d,e,n) * (n+2)/(n+1)#*k1k0(P,n+1) # . *  k1k0(dP,n) = k(dP,n+1)/k(dP,n) = (n+2)kn(P,n+1)/((n+1)kn(P,n+1) = (n+2)/(n+1)*k1k0(P,n+1)
    1/Aᴵn
end

βᴵn(P::Type{<:AbstractCOP}, n::Int) = β̃ᴵn(P, n)
function β̃ᴵn(P::Type{<:AbstractCOP}, n::Int)
    n = n - 1
    a,b,c,d,e = abcdeᴵ(P)
    Ãᴵn = Ãn(P,a,b,c,d,e,n)
    B̃ᴵn = B̃n(P,a,b,c,d,e,n)

    - B̃ᴵn / Ãᴵn
end

γᴵn(P::Type{<:AbstractCOP}, n::Int)  = γ̃ᴵn(P,n) * k1k0(P, n-1)
function γ̃ᴵn(P::Type{<:AbstractCOP}, n::Int)
    n = n - 1
    a,b,c,d,e = abcdeᴵ(P)
    Ãᴵn = Ãn(P,a,b,c,d,e,n)  # * k(dP,n+1)/k(dP,n)
    C̃ᴵn = C̃n(P,a,b,c,d,e,n)  # * k(dP,n+1)/k(dP,n-1)
    C̃ᴵn / Ãᴵn * (n+1)/n #* k1k0(P, n)     # k(dp,n)/k(dp,n-1) = (n+1)k(P,n+1)/(n k(P,n)) = (n+1)/n * k1k0(n)
end

##
##  --------------------------------------------------
##

##
## --------------------------------------------------
##
## Conversion

function Base.convert(::Type{Q},  p::P)  where {Q <: Polynomials.StandardBasisPolynomial, P <: AbstractCOP}
    p(variable(Q, p.var))
end
function Base.convert(::Type{Q},  p::P)  where {Q <: AbstractCCOP,  P <: Polynomials.StandardBasisPolynomial}
    _convert_cop(Q, p)
end

## Conversion
## * use FastTransforms, when available for T <: AbstractFloat; see connection.jl
## * use  _convert_cop when possible (needs to match σ)
## * use  convesion  through Polynomial type
function Base.convert(::Type{Q}, p::P)  where  {Q <: AbstractCCOP,  P <: AbstractCCOP}
    _convert(Q, p)
end

# work around method ambiguity introducted in abstract
# dispatch  to  specific FastTransform  method  (defined in `connection.jl`) or
# use this default
function  _convert(::Type{Q}, p::P) where {Q <: AbstractCCOP,  P <: AbstractCCOP}

    a,b,c,d,e = abcde(P)
    ā,b̄,c̄,d̄,ē = abcde(Q)

    if (a,b,c) == (ā,b̄,c̄)  &&  !(Q <: Hermite || P  <: Hermite)  # σ ==  σ̄   and not Hermite
        _convert_cop(Q, p)
    else
        T = eltype(Q)
        convert(Q, convert(Polynomial{T}, p))
    end
end

##
## --------------------------------------------------
# # scalar ops

#  avoid dispatch when N is known
function Base.:+(p::P, c::S) where {T, N, P<:AbstractCOP{T,N}, S<:Number}
    c′ =  c / k0(P) #one(T) * ⟒(P)(c)[0] #  c / k0(P), needed to add to a coefficient
    R = promote_type(promote_type(T,S), typeof(inv(k0(P))))
    iszero(c) && return (N == 0 ? zero(⟒(P){R}, p.var) :  ⟒(P){R,N}(R.(p.coeffs), p.var))
    N == 0 && return ⟒(P)(R(c), p.var)
    N == 1 && return ⟒(P)(R[p[0] + c′], p.var)
    cs = R[iszero(i) ? p[0]+c′ : p[i] for i in 0:N-1]
    return ⟒(P){R,N}(cs, p.var)
end

##
## multiply/addition/divrem with P{α…,T,N} =  Q{α...,T, M}
## We don't have (p::P{N},q::P{M}) where  {N,M, P<:AbstractCOP} as an available signature
## so we create these fall  backs, and direct +,  *, divrem to ⊕, ⊗, _divrrem  in the `register` macros
function ⊕(p::P, q::Q) where {T,N,S,M, P <: AbstractCOP{T,N}, Q <: AbstractCOP{S,M}}

    #@assert  ⟒(P) == ⟒(Q)
    #@assert eltype(p) == eltype(q)
    R = promote_type(T,S)
    Polynomials.isconstant(p)  && return q + p(0)
    Polynomials.isconstant(q)  && return p + q(0)

    p.var != q.var && throw(ArgumentError("Variables don't  match"))

    if N==M
        as = [p[i]+q[i] for i in 0:N-1]
        return ⟒(P)(as, p.var) # might have trailing zeros here
    end

    N > M ? ⊕′(p,q) : ⊕′(q,p)

end

# assumes N >  M
@generated  function ⊕′(p1::P, p2::Q) where {T,N,S,M, P<:AbstractCOP{T,N}, Q<:AbstractCOP{S,M}}

    #@assert  ⟒(P) == ⟒(Q)
    R = promote_type(T,S)

    exprs = Any[nothing for i = 1:N]
    for i in  1:M
        exprs[i] = :(p1.coeffs[$i] + p2.coeffs[$i])
    end
    for i in M+1:N
        exprs[i] =:(p1.coeffs[$i])
    end

    return quote
        Base.@_inline_meta
        ⟒(P){$(R),$(N)}(tuple($(exprs...)), p1.var)
    end

end

function ⊗(p::P, q::Q) where {P <: AbstractCOP, Q <: AbstractCOP}

    Polynomials.isconstant(p)  && return q * p(0)
    Polynomials.isconstant(q)  && return p * q(0)
    p.var != q.var && throw(ArgumentError("Variables don't  match"))

    # use connection for linearization;  note:  evalauation  is  faster than _convert_cop
    p′,q′ = _convert_cop.(Polynomial, (p,q))
    _convert_cop(⟒(P), p′ * q′)
#    convert(⟒(P), convert(Polynomial, p) * convert(Polynomial, q))

end

## Modifications needed due to `N` in the type parameter

function Polynomials.truncate(p::P;
                               rtol::Real = Base.rtoldefault(real(T)),
                               atol::Real = 0,) where {T,N,P<:AbstractCOP{T,N}}
    ps = coeffs(p)
    max_coeff = maximum(abs, ps)
    thresh = max_coeff * rtol + atol
    map!(c->abs(c) <= thresh ? zero(T) : c, ps,ps)
    ⟒(P){T}(ps, p.var)
end

Polynomials.truncate!(p::P;
                      rtol::Real = Base.rtoldefault(real(T)),
                      atol::Real = 0,) where {T,N,P<:AbstractCOP{T,N}} = error("`truncate!` not defined")

function Base.chop(p::P;
                   rtol::Real = Base.rtoldefault(real(T)),
                   atol::Real = 0,) where {T,N,P<:AbstractCOP{T,N}}

    N == 0 && return p
    i = N-1
    while i >= 0
        val = p[i]
        if !isapprox(val, zero(T); rtol = rtol, atol = atol)
            break
        end
        i -= 1
    end
    ⟒(P)(coeffs(p)[1:i+1], p.var)
end

Polynomials.chop!(p::P;
                  rtol::Real = Base.rtoldefault(real(T)),
                  atol::Real = 0,) where {T,N,P<:AbstractCOP{T,N}} = error("`chop!` not defined")

# use pn= [â,b̂,ĉ] ⋅ [p'_{n+1}, p'_n, p'_{n-1}] to
# find expression for p' in terms of p
function Polynomials.derivative(p::P, order::Integer=1) where {P <:AbstractCOP}

    R = eltype(one(eltype(p))/1)
    d = degree(p)
    order < 0 && throw(ArgumentError("order must  be non-negative"))
    order == 0 && return ⟒(P){R}(coeffs(p), p.var)
    d < order && return zero(⟒(P){R})

    as = zeros(R, d)
    ps = R.(coeffs(p))
    for n = d-1:-1:1
        a,b,c = ân(P,n),b̂n(P,n),ĉn(P,n)
        if !iszero(a)
            pn = ps[1+n+1]
            as[1+n] = pn/a
            ps[1+n] -= pn*b/a
            ps[1+n-1] -= pn*c/a
        end
    end
    a,b = ân(P,0),b̂n(P,0)
    p1 = ps[1+1]
    as[1+0] = p1/a

    dp = ⟒(P)(as, p.var)
    order == 1 ? dp : derivative(dp, order - 1)

end

function Polynomials.integrate(p::P, C::Number=0) where {P <: AbstractCOP}

    T,S = eltype(p), typeof(C)
    R = promote_type(typeof(one(T) / 1), S)
    Q = ⟒(P){R}
    if hasnan(p) || isnan(C)
        return   Q(NaN)
    end

    n = degree(p)
    if n == -1
        return zero(Q, p.var)
    elseif n == 0
        return C*one(Q, p.var)  + p(0)*variable(Q, p.var)
    end

    as = zeros(R, n + 2)

    # case d=0 we do by hand,  as
    # P_0(x) = c, so ∫P_o = x = c*variable(P)
    c₀,c₁ = coeffs(variable(p))
    pd = first(p.coeffs)
    as[1] = pd * c₀
    as[2] = pd * (c₁ * k0(Q))
    @inbounds for d in 1:n
        pd = p.coeffs[d+1]
        as[1 + d + 1] += pd * ân(Q, d)
        as[1 + d]     += pd * b̂n(Q, d)
        if  d > 1
            as[1 + d - 1] += pd * ĉn(Q, d)
        end
    end

    # adjust constant
    ∫p = Q(as,  p.var)
    return ∫p - ∫p(0) + Q(C)

end


text:
# greedy approach to fit
struct GreedyIntegrator_Cache
    Y
    X
    XZ
    XV
    XU
    ZIntegrator
end

struct GreedyIntegrator_Params
    Z_alg
    Z_kwargs
end

struct GreedyIntegrator <: AbstractDLRAlgorithm
    alg_params::GreedyIntegrator_Params
end
function GreedyIntegrator(; Z_alg = Tsit5(), Z_kwargs = Dict{Symbol,Any}())
    params = GreedyIntegrator_Params(Z_alg, Z_kwargs)
    return GreedyIntegrator(params)
end

function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::TwoFactorRepresentation, dt)
    X = zeros(size(u))
    r = rank(u)
    n = size(X,1)
    XZ = zeros(n,r)
    return GreedyIntegrator_Cache(prob.y, X, XZ, nothing, nothing, nothing)
end

function alg_cache(prob::MatrixDataProblem, alg::GreedyIntegrator, u::SVDLikeRepresentation, dt)
    X = zeros(size(u))
    r = rank(u)
    n, m = size(X)
    XV = zeros(n,r)
    XU = zeros(m,r)
    return GreedyIntegrator_Cache(prob.y, X, nothing, XV, XU, nothing)
end

function alg_cache(prob::MatrixHybridProblem, alg::GreedyIntegrator, u, dt)
    XZ = zeros(size(u,1), rank(u))
    X = zeros(size(u)...)
    ZProblem = ODEProblem(prob.f, u.Z, prob.tspan, u.U)
    ZIntegrator = init(ZProblem, alg.alg_params.Z_alg; save_everystep=false, alg.alg_params.Z_kwargs...)
    return GreedyIntegrator_Cache(prob.y, X, XZ, nothing, nothing, ZIntegrator)
end

function init(prob::MatrixDataProblem, alg::GreedyIntegrator, dt)
    t0, tf = prob.tspan
    @assert tf > t0 "Integration in reverse time direction is not supported"
    u = deepcopy(prob.u0)
    # initialize solution
    sol = init_sol(dt, t0, tf, prob.u0)
    # initialize cache
    cache = alg_cache(prob, alg, u, dt)
    sol.Y[1] = deepcopy(prob.u0)
    return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0)
end

function init(prob::MatrixHybridProblem, alg::GreedyIntegrator, dt)
    t0, tf = prob.tspan
    @assert tf > t0 "Integration in reverse time direction is not supported"
    u = deepcopy(prob.u0)
    # initialize solution
    sol = init_sol(dt, t0, tf, prob.u0)
    # initialize cache
    cache = alg_cache(prob, alg, u, dt)
    sol.Y[1] = deepcopy(prob.u0)
    return DLRIntegrator(u, t0, dt, sol, alg, cache, typeof(prob), 0)
end

function greedy_step!(u::TwoFactorRepresentation, cache, t, dt, ::Type{<:MatrixHybridProblem})
    @unpack Z, U = u
    @unpack X, Y, ZIntegrator, XZ = cache
    step!(ZIntegrator, dt, true)
    Z .= ZIntegrator.u
    update_data!(X, Y, t, dt)
    mul!(XZ, X, Z)
    Q, _, P = svd(XZ)
    mul!(U, Q, P')
end

function greedy_step!(u::TwoFactorRepresentation, cache, t, dt, ::Type{<:MatrixDataProblem})
    @unpack Z, U = u
    @unpack Y, X, XZ = cache
    update_data!(X,Y,t,dt)
    mul!(Z, X', U)
    mul!(XZ, X, Z)
    Q, _, P = svd(XZ)
    mul!(U, Q, P')
end

function greedy_step!(u::SVDLikeRepresentation, cache, t, dt, ::Type{<:MatrixDataProblem})
    @unpack U, S, V = u
    @unpack Y, X, XV, XU = cache
    update_data!(X,Y,t,dt)
    mul!(XV,X,V)
    mul!(XU,X',U)
    U .= Matrix(qr(XV).Q)
    V .= Matrix(qr(XU).Q)
    mul!(XV,X,V)
    mul!(S,U',XV)
end

function step!(integrator::DLRIntegrator, ::GreedyIntegrator, dt)
    @unpack u, t, iter, cache, probType = integrator
    greedy_step!(u, cache, t, dt, probType)
    integrator.t += dt
    integrator.iter += 1
end


text:
using Documenter, SineWaves

makedocs(modules = [SineWaves],
         sitename = "SineWaves.jl",
         pages = Any[
            "Home" => "index.md",
            "Wizard" => "wizard.md",
            "BinaryBuilder" => "builder.md",
            "C calls" => "ccalls.md",
            "Functions" => "api.md",
            "Examples" => "examples.md"
        ])

deploydocs(
    repo = "github.com/jakubwro/SineWaves.jl.git",
    target = "build"
)


text:
"""
This function returns a vector of investors in the simulation.
"""
function create_investors(simulation_data::AgentSimulationData)
    horizon = get_total_horizon(get_case(simulation_data))
    simulation_data_dir = get_data_dir(get_case(simulation_data))
    dir_name = joinpath(get_data_dir(get_case(simulation_data)), "investors")
    investor_names = readdir(dir_name)
    investors = Vector{Investor}(undef, length(investor_names))
    for i = 1:length(investor_names)
        investor_dir = joinpath(dir_name, "$(investor_names[i])")

        rep_hour_weight = get_rep_hour_weight(simulation_data)

        # Read Investor Characteristics
        characteristics = read_data(joinpath(investor_dir, "characteristics.csv"))

        forecast_type = get_forecast_type(get_case(simulation_data))

        if forecast_type == "perfect"
            parameter_values = [get_annual_growth(simulation_data) for i in 1:horizon]
            scenario_data = [Scenario("perfect", 1.0, nothing, parameter_values)]
            forecast = Perfect(scenario_data)

        elseif forecast_type == "imperfect"

            # Read belief data file
            if get_info_symmetry(get_case(simulation_data))
                belief_filename = joinpath(simulation_data_dir, "markets_data", "symmetric_belief.csv")
            else
                belief_filename = joinpath(investor_dir, "markets_data", "investor_belief.csv")
            end

            @assert isfile(belief_filename)

            investor_belief_data = read_data(belief_filename)

            param_names= investor_belief_data.parameters

            @assert in("initial_estimate", names(investor_belief_data))
            initial_estimate = AxisArrays.AxisArray(investor_belief_data.initial_estimate, param_names)   # initial state initial_estimate

            # Create Kalman Filter for Belief Update
            if get_belief_update(get_case(simulation_data))
                @assert in("initial_error_cov", names(investor_belief_data))
                @assert in("process_cov", names(investor_belief_data))
                @assert in("measurement_cov", names(investor_belief_data))

                initial_error_covariance = AxisArrays.AxisArray(collect(LinearAlgebra.Diagonal(investor_belief_data.initial_error_cov)), param_names, param_names) # initial error covariance

                Q = AxisArrays.AxisArray(investor_belief_data.process_cov, param_names)        # process covariance
                R = AxisArrays.AxisArray(investor_belief_data.measurement_cov, param_names)    # measurement covariance

                belief = InvestorBelief(Q, R)
                kf = KalmanFilter(belief, initial_estimate, initial_error_covariance)

            else
                kf = nothing
            end

            # Populate scenario data
            scenario_data = Scenario[]
                # If user has provided parameter values for each scenario

            if get_uncertainty(get_case(simulation_data))
                # Assert that user has provided parameter multiplier values for each scenario
                if get_info_symmetry(get_case(simulation_data))
                    scenario_file_name = joinpath(simulation_data_dir, "markets_data", "symmetric_scenario_multiplier_data.csv")
                else
                    scenario_file_name = joinpath(investor_dir, "markets_data", "scenario_multiplier_data.csv")
                end

                @assert isfile(scenario_file_name)
                scenario_df = read_data(scenario_file_name)
                num_scenarios = DataFrames.nrow(scenario_df)
                for s in 1:num_scenarios
                    name = scenario_df[s, :scenario]
                    probability = scenario_df[s, :probability]
                    parameter_multipliers = Dict{String, Float64}()
                    parameter_values = [AxisArrays.AxisArray(zeros(length(param_names), horizon), param_names, collect(1:horizon)) for i in 1:horizon]
                    for param in param_names
                        if in(param, names(scenario_df))
                            parameter_multipliers[param] = scenario_df[s, Symbol(param)]
                            for i in 1:horizon
                                parameter_values[1][param, i] = parameter_multipliers[param] * initial_estimate[param]
                            end
                        else
                            parameter_multipliers[param] = 1.0
                            for i in 1:horizon
                                parameter_values[1][param, i] = initial_estimate[param]
                            end
                        end
                    end
                    push!(scenario_data, Scenario(name, probability, parameter_multipliers, parameter_values))
                end
            else
                name = "scenario_1"
                probability = 1.0
                parameter_multipliers = Dict(param => 1.0 for param in param_names)
                parameter_values = [AxisArrays.AxisArray(zeros(length(param_names), horizon), param_names, collect(1:horizon)) for i in 1:horizon]
                for param in param_names
                    for i in 1:horizon
                        parameter_values[1][param, i] = initial_estimate[param]
                    end
                end

                push!(scenario_data, Scenario(name, probability, parameter_multipliers, parameter_values))
            end

            forecast = Imperfect(kf, scenario_data)

        end
        #Empty vector of projects.
        projects = Project{<:BuildPhase}[]

        projectdata_existing = extract_projectdata(investor_dir, "projectexisting.csv")
        projectdata_options = extract_projectdata(investor_dir, "projectoptions.csv")

        sys_UC = get_system_UC(simulation_data)

        #Append existing and option projects.
        append!(projects, create_project_existing(projectdata_existing,
                                                  simulation_data,
                                                  sys_UC,
                                                  investor_names[i],
                                                  investor_dir,
                                                  get_name.(scenario_data)))

        append!(projects, create_project_options(projectdata_options,
                                                      simulation_data,
                                                      investor_names[i],
                                                      investor_dir,
                                                      get_name.(scenario_data)))

        add_investor_project_availability!(simulation_data_dir, projects)

        #Names of the markets in which the investor is participating.
        markets = Symbol[]
        simulation_markets = get_markets(simulation_data)

        for m in keys(simulation_markets)
            if simulation_markets[m]
                if m == :Reserves
                    push!(markets, :ReserveUp)
                    push!(markets, :ReserveDown)
                else
                    push!(markets, m)
                end
            end
        end

        # Empty market prices struct
        market_prices = MarketPrices()

        capital_cost_multiplier = characteristics.capital_cost_multiplier[1]
        max_annual_projects = characteristics.max_annual_projects[1]

        # Risk Preference
        if get_risk_aversion(get_case(simulation_data))
            risk_preference_type = lowercase(characteristics.risk_preference[1])
            if risk_preference_type == "neutral"
                risk_preference = RiskNeutral()
            elseif risk_preference_type == "averse"
                risk_preference = RiskAverse(characteristics.uf_constant[1], characteristics.uf_multiplier[1], characteristics.uf_risk_coefficient[1])
            end
        else
            risk_preference = RiskNeutral()
        end

        investors[i] = Investor(investor_names[i],
                                investor_dir,
                                projects,
                                markets,
                                market_prices,
                                rep_hour_weight,
                                forecast,
                                capital_cost_multiplier,
                                max_annual_projects,
                                risk_preference)
    end

    return investors
end