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

using Markdown
using InteractiveUtils

# ╔═╡ 96265874-5d54-11eb-1eb0-1f1f4f0b55d2
begin
	using Pkg
	Pkg.add("DataFrames")
	Pkg.add("CSV")
	Pkg.add("Plots")
end

# ╔═╡ d105ceb6-5d54-11eb-13b4-63cbff4d3770
begin
	using Dates
	using CSV
	using DataFrames
	using Plots
end

# ╔═╡ c76a5d2e-639c-11eb-0184-edbef6db4ef2
md"*Выполнили:* Скворцов Иван, Никулина Евгения, Кордзахия Натела, Татаринов Артем (БЭК-181)"

# ╔═╡ 4fe84764-5d5e-11eb-3024-31e0b9505741
md"# Используемые пакеты"

# ╔═╡ aabae65e-5d60-11eb-1bf6-1b1db9ee1463
md"# Вводные слова

В данном обзоре функционала Julia мы будем использовать две таблицы, содержащие обезличенные данные одной известной сети магазинов одежды.

Таблица `transactions_footfall` содержит данные по количеству посетителей и транзакций в магазинах сети.

Таблица `events` содержит количество кликов (событий) в приложении, которое продавцы-консультанты используют для работы с клиентами.

Целью нашей работы станет исследование зависимости между **частотой использования приложения** и **конверсией** (отношением кол-ва транзакций к общему числу клиентов).
"

# ╔═╡ 61daeccc-5d5e-11eb-15ff-dd2ea6cace20
md"# Подгрузка датасетов

Данные хранятся в формате `csv`. Для их корректной подгрузки мы задаем формат дат (переменная `dftm`), а также тип первого столбца -- `String`, поскольку ID магазина является категориальной переменной.
"

# ╔═╡ 2a6b4728-5d60-11eb-2b52-fd202a979ad0
dfmt = dateformat"yyyy-mm-dd"

# ╔═╡ d0ed690c-5d54-11eb-317f-292766404960
transactions_footfall = CSV.read("transactions_footfall.csv", DataFrame, types=Dict(1=>String), dateformat=dfmt) # объект Dict в аргументе types задает желательные форматы столбцов

# ╔═╡ d0d19420-5d54-11eb-3a88-b72c160f299e
usage = CSV.read("usage.csv", DataFrame, types=Dict(1=>String), dateformat=dfmt)

# ╔═╡ 7b9ab988-5d62-11eb-240e-d5e0188871ff
md"
# Первичное исследование данных

При помощи функции `describe` выясняем, что имеющиеся данные захватывают период с 1 по 30 ноября 2020 года."

# ╔═╡ 73789694-5d62-11eb-39f3-1b71b1387a03
describe(transactions_footfall)

# ╔═╡ a85b4b34-5d5f-11eb-3875-aba54752e12b
describe(usage)

# ╔═╡ 783d8f8a-5d63-11eb-30a3-fba836f760c5
md"При этом данные по пользованию приложением есть не для всех магазинов: скорее всего, в некоторых из них оно не использовалось вовсе. Это дает нам возможность посмотреть на работу различных типов объединения данных: *left join*, *right join* и *inner join*."

# ╔═╡ fb527712-5d62-11eb-0b18-c57c96010935
length(unique(transactions_footfall.store_id))

# ╔═╡ 2249e9ac-5d63-11eb-1e15-adbc290ab7df
length(unique(usage.store_id))

# ╔═╡ b3e610ec-5d64-11eb-0aaa-fbb9619a4936
md"При помощи операций над множествами мы можем выяснить, насколько полными будут данные после их объединения. Для этого зададим множества значений столбца `store_id` в для обоих датасетов:"

# ╔═╡ 42dc4646-5d64-11eb-2baf-015f759a2996
begin
	stores_in_tf = Set(transactions_footfall.store_id)
	stores_in_ug = Set(usage.store_id)
end

# ╔═╡ 1af80bf0-5d65-11eb-0d17-63353d720eb1
md"В частности, для 76 из 77 магазинов, использовавших приложение в ноябре, имеются данные по трафику и транзакциям:"

# ╔═╡ 3e4dd1ee-5d64-11eb-021e-8f34c40da0a3
intersect(stores_in_tf, stores_in_ug)

# ╔═╡ 653867be-5d65-11eb-2091-89e3ce947048
md"17 магазинов не использовали приложение, но обслуживали посетителей:"

# ╔═╡ 66165a06-5d65-11eb-14b1-d996eb60e56c
setdiff(stores_in_tf, stores_in_ug)

# ╔═╡ b5c9dc8a-5d65-11eb-1e5f-71e4771e5409
md"Опять же, у 1 магазина из использовавших приложение нет данных по числу посетителей (возможно, он не работал, а приложение использовалось для внутренних процессов)."

# ╔═╡ 9b592126-5d65-11eb-27bd-adbaf0581c61
setdiff(stores_in_ug, stores_in_tf)

# ╔═╡ 52080d92-5d66-11eb-245f-91638f919f6d
md"Можем вывести данные по пользованию для этого магазина. Видим, что приложение использовалось всего четыре дня."

# ╔═╡ f2e872b6-5d65-11eb-1913-4503dcea7a9e
usage[usage[!, "store_id"] .== "4659", :]

# ╔═╡ f7076286-5d69-11eb-3005-e9cdc1e55211
md"Наконец, в целом две имеющиеся таблицы покрывают 94 магазина сети:"

# ╔═╡ 078d3798-5d6a-11eb-09d9-87da93fda83f
union(stores_in_ug, stores_in_tf)

# ╔═╡ bd483046-5d66-11eb-22b5-b9b0386b9692
md"# Объединение таблиц (join)

Объединение таблиц обычно понимается в смысле join'ов в том виде, в каком они представлены в языке SQL. Ниже представлено схематическое представление четырех основных типов объединений. В этом разделе мы разберем каждый из них на базе функционала Julia DataFrames.

![joins](https://orkhanalyshov.com/media/uploads/files/joins.png)
"

# ╔═╡ 9b10661a-5d6b-11eb-2695-b32e5d8e70e4
md"## Синтаксис

Объединение таблиц в Julia производится при помощи набора функций, полный перечень которых представлен в [документации](https://dataframes.juliadata.org/stable/man/joins/). Разберем общую логику их синтаксиса.

Функции вида `*join(df1, df2, on = ...)` принимают три обязательных аргумента: *левая (первая) таблица*, *правая (вторая) таблица* и *столбцы-индексы (key columns)*.

Основная сложность может возникнуть при задании столбцов-индексов. Разберем основные случаи, которые встречаются при работе с реальными данными:
- Названия столбцов-индексов **совпадают** в обеих таблицах.
  - Один индекс: `on = :idx_col_name`
  - Несколько индексов: `on = [:idx_col_name1, :idx_col_name2, ...]`
- Названия столбцов-индексов **не совпадают** в обеих таблицах. В таком случае нужно задать соответствие между названиями столбцов в датасетах:
  - Первый способ, кортежи: `on = [(:idx_1_df1, :idx_1_df2), (:idx_2_df1, :idx_2_df2), ...]`
  - Второй способ, словари: `on = [:idx_1_df1 => :idx_1_df2, :idx_2_df1 => :idx_2_df2, ...]`

Иногда бывает необходимо проверить, являются ли сочетания индексов уникальными в каждом из датасетов (например, убедиться, что в таблице встречается лишь одно значение  переменной для пары store_id-date). В таком случае используется дополнительный аргумент функции `*join` -- `validate`, принимающий кортеж из двух логических значений true/false. Первое значение отвечает за проверку уникальности в левом датасете, второе значение -- за проверку уникальности в правом датасете. Если будут обнаружены повторяющиеся индексы, функция выдаст ошибку типа

`ERROR: ArgumentError: Merge key(s) are not unique in both df1 and df2. First duplicate in df1 at 3. First duplicate in df2 at 3`
"

# ╔═╡ ed38d744-5d67-11eb-2c2b-17b39daad6bd
md"## Inner join

Этот тип объединения подразумевает, что в итоговый датасет включаются только те наблюдения, индексы которых содержатся в обеих таблицах. В нашем случае мы имеем дело с двойным индексом: нам необходимо объединять данные по столбцам *store_id* и *date* одновременно."

# ╔═╡ ecdc3856-5d67-11eb-199e-f97adc1bfcee
inner = innerjoin(usage, transactions_footfall, on = [:store_id, :date])

# ╔═╡ d08869f8-5d68-11eb-2be4-a374d299b641
md"Исходя из нашего знания данных, полученного в предыдущем разделе, мы ожидаем, что inner join позволит получить объединенный датасет для **76 магазинов**, у которых имеются данные и по пользованию, и по посетителям. Проверим, что это действительно так:"

# ╔═╡ ec93af50-5d67-11eb-05a0-1db28fbc0c8d
Set(inner.store_id)

# ╔═╡ ec0af2fa-5d67-11eb-32bc-214e41ec47a8
md"## Full join

Этот тип объединения сохранит все индексы, которые имеются в двух таблицах. При этом недостающие данные будут заполнены пустыми значениями. В Julia соответствующая функция называется `outerjoin` (не путать с настоящим outer join, сохраняющим только наблюдения, индексы которых встречаются лишь в одном из датасетов!)"

# ╔═╡ 44fa01de-5d69-11eb-2c65-9d8caeb1112a
full = outerjoin(usage, transactions_footfall, on = [:store_id, :date])

# ╔═╡ 44e45802-5d69-11eb-3744-b7da7627a568
md"Как и следовало ожидать, полученный датасет включает данные по 94 магазинам:"

# ╔═╡ 44cff20e-5d69-11eb-200e-c1e5ce7c7f84
Set(full.store_id)

# ╔═╡ 44bbc52c-5d69-11eb-31bd-d37a4ea8499d
md"## Left join

Этот тип объединения сохранит все индексы, которые содержатся в левой (первой) таблице. Недостающие данные из правой таблицы будут заполнены пустыми значениями."

# ╔═╡ 44693b72-5d69-11eb-3e95-43a9f012adcb
left = leftjoin(usage, transactions_footfall, on = [:store_id, :date])

# ╔═╡ 8455b2d2-5d6a-11eb-3749-7f09bc737d87
md"Разумеется, полученная таблица содержит данные по 77 магазинам:"

# ╔═╡ 84434554-5d6a-11eb-3f58-71f79d52237e
Set(left.store_id)

# ╔═╡ 842d68d6-5d6a-11eb-10a4-0d2321571d6e
md"## Right join

Этот тип объединения сохранит все индексы, которые содержатся в правой (второй) таблице. Недостающие данные из левой таблицы будут заполнены пустыми значениями."

# ╔═╡ 83e1019c-5d6a-11eb-38b0-e19f9f271e03
right = rightjoin(usage, transactions_footfall, on = [:store_id, :date])

# ╔═╡ bf10e93c-5d6a-11eb-277e-1de238f1c555
md"Полученная таблица содержит данные по 93 магазинам:"

# ╔═╡ a89c8efe-5d6a-11eb-2860-6767e9b76031
Set(right.store_id)

# ╔═╡ a86f1c5a-5d6a-11eb-13b5-47f7b59ea614
md"Поскольку left и right объединения относительны, right join также можно сделать, поменяв порядок таблиц внутри функции `leftjoin`. "

# ╔═╡ a857f938-5d6a-11eb-2130-73d6529537d9
right_alternative = leftjoin(transactions_footfall, usage, on = [:store_id, :date])

# ╔═╡ a8412eb0-5d6a-11eb-1842-ed9a2583408c
Set(right_alternative.store_id)

# ╔═╡ 79301162-5d6b-11eb-16b2-a19850841803
md"Более подробную информацию об объединениях таблиц в Julia DataFrames можно найти [здесь](https://dataframes.juliadata.org/stable/man/joins/)."

# ╔═╡ 2583d2b0-5d6b-11eb-35e1-35b483085f3e
md"# Конкатенация таблиц

Конкатенация таблицы - это простое склеивание вдоль одной из осей. Ниже мы рассмотрим конкатенацию по горизонтали - `hcat()` и по вертикали - `vcat()`. Обе функции происходят от базовой функции `cat(..., dims = (m,n))`, которую удобно использовать, если надо объединить таблицы с тремя и более измерениями. Для тех, кто знаком с питоном - эта фунция полностью аналогична `pd.concat()`

![joins](https://i.ibb.co/NFSP4qf/image.png)"

# ╔═╡ 350a83b2-5d70-11eb-3d4a-cd801bc268ec
md"## По вертикали

Конкатенация по вертикали позволяет объединить два датасета, если у них совпадает количество столбцов (в нашем примере три), при этом количество строк для объединения может быть любым:"

# ╔═╡ 47768ae0-5d73-11eb-1bff-a168c868ef07
a = usage[1:4, :]

# ╔═╡ 46bbd506-5d73-11eb-03c9-73316a55f932
b = usage[5:8, :]

# ╔═╡ 34f88772-5d70-11eb-0948-e19ff61d0231
v_conc = vcat(a,b)

# ╔═╡ 34dfb7b0-5d70-11eb-2807-31dde3f1f3d2
md"## По горизонтали

Конкатенация по горизонтали позволяет объединить два датасета, если у них совпадает количество строк. При этом количество столбцов для объединения может быть любым:"

# ╔═╡ 4622f334-5d73-11eb-2277-cd5a6049d0db
c = usage[:, 1:1]

# ╔═╡ 4611174c-5d73-11eb-1db8-65f15c3841a3
d = usage[:, 2:3]

# ╔═╡ 45fa7816-5d73-11eb-2f77-3b0e0e10fa5f
h_conc = hcat(c,d)

# ╔═╡ f5ab80be-5fbc-11eb-3631-05292cb739c5
md"В [документации](https://docs.julialang.org/en/v1/base/arrays/#Concatenation-and-permutation) можно подробнее прочитать о конкатенации в Julia."

# ╔═╡ 348a414a-5d70-11eb-3ae9-fffe5497977e
md"# Корреляционный анализ

Применяя полученные знания об объединении таблиц, выясним характер взаимосвязи между **частотой использования приложения** и **конверсией**.

Поскольку для этого нам нужны наиболее полные данные, будем использовать inner join для создания рабочего датасета:
"

# ╔═╡ cfe2cc5a-5d54-11eb-0d91-b568ed9aae7e
data = innerjoin(usage, transactions_footfall, on = [:store_id, :date])

# ╔═╡ 916d88e2-5d70-11eb-2b1f-d529130b1155
md"Исключим наблюдения с пропущенными значениями (они встречаются только в столбце `footfall`):"

# ╔═╡ a3dae6e4-5d57-11eb-34b9-b723c2639ce5
dropmissing!(data)

# ╔═╡ 364ac9cc-5d71-11eb-0047-8595df7520f3
md"Исключим наблюдения с нулевым количеством посетителей или транзакций: скорее всего, в такие дни магазин не работал."

# ╔═╡ 35f773bc-5d71-11eb-169f-6fc9ff508c55
begin
	filter!(row -> row.footfall != 0, data)
	filter!(row -> row.transactions != 0, data)
end

# ╔═╡ aa14d1e6-5d70-11eb-2d86-ed9f13e33513
md"Зададим новые переменные: `conversion` (конверсия) и `CPP` (clicks per person, число кликов в приложении на посетителя)."

# ╔═╡ 86ddd0a0-5d58-11eb-38ff-bb9980ad61f6
begin
	data[!, "conversion"] = data.transactions ./ data.footfall
	data[!, "CPP"] = data.events ./ data.footfall
end

# ╔═╡ a3efaaca-5d57-11eb-3d23-c535e3b971d6
data

# ╔═╡ d6d2d6ce-5d70-11eb-3e56-1d24a0f7d104
md"Отфильтруем наблюдения, руководствуясь тем, что конверсия не может превышать 1, а число кликов на посетителя вряд ли превысит 25."

# ╔═╡ 6faf1d30-5d5d-11eb-2bee-75b8bdf318c1
data_clean = data[(data[!, "conversion"] .< 1) .& (data[!, "CPP"] .< 10), :]

# ╔═╡ 0982d144-5d73-11eb-14a5-cf57c0fb38d1
plot(data_clean.CPP, data_clean.conversion, seriestype = :scatter, title = "Clicks-per-person vs. Conversion", dpi = 400, yformatter = y -> string(Int64(floor(y * 100)), "%"), smooth = true, label = "", xlabel = "Clicks-per-person", ylabel = "Conversion")

# ╔═╡ 25bbd8fe-639c-11eb-2507-9dac2bb2e31a
md"Можно заметить, что наблюдается положительная зависимость между интенсивностью пользования приложением и конверсией! Круто!"

# ╔═╡ 67b97322-5d72-11eb-0911-fbfc91084b23
md"**Гетероскедастичность -- это отдельный разговор :)**"

# ╔═╡ Cell order:
# ╟─c76a5d2e-639c-11eb-0184-edbef6db4ef2
# ╟─4fe84764-5d5e-11eb-3024-31e0b9505741
# ╠═96265874-5d54-11eb-1eb0-1f1f4f0b55d2
# ╠═d105ceb6-5d54-11eb-13b4-63cbff4d3770
# ╟─aabae65e-5d60-11eb-1bf6-1b1db9ee1463
# ╟─61daeccc-5d5e-11eb-15ff-dd2ea6cace20
# ╟─2a6b4728-5d60-11eb-2b52-fd202a979ad0
# ╠═d0ed690c-5d54-11eb-317f-292766404960
# ╠═d0d19420-5d54-11eb-3a88-b72c160f299e
# ╟─7b9ab988-5d62-11eb-240e-d5e0188871ff
# ╠═73789694-5d62-11eb-39f3-1b71b1387a03
# ╠═a85b4b34-5d5f-11eb-3875-aba54752e12b
# ╟─783d8f8a-5d63-11eb-30a3-fba836f760c5
# ╠═fb527712-5d62-11eb-0b18-c57c96010935
# ╠═2249e9ac-5d63-11eb-1e15-adbc290ab7df
# ╟─b3e610ec-5d64-11eb-0aaa-fbb9619a4936
# ╠═42dc4646-5d64-11eb-2baf-015f759a2996
# ╟─1af80bf0-5d65-11eb-0d17-63353d720eb1
# ╠═3e4dd1ee-5d64-11eb-021e-8f34c40da0a3
# ╟─653867be-5d65-11eb-2091-89e3ce947048
# ╠═66165a06-5d65-11eb-14b1-d996eb60e56c
# ╟─b5c9dc8a-5d65-11eb-1e5f-71e4771e5409
# ╠═9b592126-5d65-11eb-27bd-adbaf0581c61
# ╟─52080d92-5d66-11eb-245f-91638f919f6d
# ╠═f2e872b6-5d65-11eb-1913-4503dcea7a9e
# ╟─f7076286-5d69-11eb-3005-e9cdc1e55211
# ╠═078d3798-5d6a-11eb-09d9-87da93fda83f
# ╟─bd483046-5d66-11eb-22b5-b9b0386b9692
# ╟─9b10661a-5d6b-11eb-2695-b32e5d8e70e4
# ╟─ed38d744-5d67-11eb-2c2b-17b39daad6bd
# ╠═ecdc3856-5d67-11eb-199e-f97adc1bfcee
# ╟─d08869f8-5d68-11eb-2be4-a374d299b641
# ╠═ec93af50-5d67-11eb-05a0-1db28fbc0c8d
# ╟─ec0af2fa-5d67-11eb-32bc-214e41ec47a8
# ╠═44fa01de-5d69-11eb-2c65-9d8caeb1112a
# ╟─44e45802-5d69-11eb-3744-b7da7627a568
# ╠═44cff20e-5d69-11eb-200e-c1e5ce7c7f84
# ╟─44bbc52c-5d69-11eb-31bd-d37a4ea8499d
# ╠═44693b72-5d69-11eb-3e95-43a9f012adcb
# ╟─8455b2d2-5d6a-11eb-3749-7f09bc737d87
# ╠═84434554-5d6a-11eb-3f58-71f79d52237e
# ╟─842d68d6-5d6a-11eb-10a4-0d2321571d6e
# ╠═83e1019c-5d6a-11eb-38b0-e19f9f271e03
# ╟─bf10e93c-5d6a-11eb-277e-1de238f1c555
# ╠═a89c8efe-5d6a-11eb-2860-6767e9b76031
# ╟─a86f1c5a-5d6a-11eb-13b5-47f7b59ea614
# ╠═a857f938-5d6a-11eb-2130-73d6529537d9
# ╠═a8412eb0-5d6a-11eb-1842-ed9a2583408c
# ╟─79301162-5d6b-11eb-16b2-a19850841803
# ╟─2583d2b0-5d6b-11eb-35e1-35b483085f3e
# ╟─350a83b2-5d70-11eb-3d4a-cd801bc268ec
# ╠═47768ae0-5d73-11eb-1bff-a168c868ef07
# ╠═46bbd506-5d73-11eb-03c9-73316a55f932
# ╠═34f88772-5d70-11eb-0948-e19ff61d0231
# ╟─34dfb7b0-5d70-11eb-2807-31dde3f1f3d2
# ╠═4622f334-5d73-11eb-2277-cd5a6049d0db
# ╠═4611174c-5d73-11eb-1db8-65f15c3841a3
# ╠═45fa7816-5d73-11eb-2f77-3b0e0e10fa5f
# ╟─f5ab80be-5fbc-11eb-3631-05292cb739c5
# ╟─348a414a-5d70-11eb-3ae9-fffe5497977e
# ╠═cfe2cc5a-5d54-11eb-0d91-b568ed9aae7e
# ╟─916d88e2-5d70-11eb-2b1f-d529130b1155
# ╠═a3dae6e4-5d57-11eb-34b9-b723c2639ce5
# ╟─364ac9cc-5d71-11eb-0047-8595df7520f3
# ╠═35f773bc-5d71-11eb-169f-6fc9ff508c55
# ╟─aa14d1e6-5d70-11eb-2d86-ed9f13e33513
# ╠═86ddd0a0-5d58-11eb-38ff-bb9980ad61f6
# ╠═a3efaaca-5d57-11eb-3d23-c535e3b971d6
# ╟─d6d2d6ce-5d70-11eb-3e56-1d24a0f7d104
# ╠═6faf1d30-5d5d-11eb-2bee-75b8bdf318c1
# ╠═0982d144-5d73-11eb-14a5-cf57c0fb38d1
# ╟─25bbd8fe-639c-11eb-2507-9dac2bb2e31a
# ╟─67b97322-5d72-11eb-0911-fbfc91084b23

max_stars_repo_path: Table joins/Table joins.jl

max_issues_repo_path: Table joins/Table joins.jl

max_forks_repo_path: Table joins/Table joins.jl


text:
function simplify_constants(O::Operation)
    O_last = nothing
    _O = O
    while _O != O_last
        O_last = _O
        _O = _simplify_constants(_O)
        if typeof(_O) <: Operation
            _O = Operation(_O.op,simplify_constants.(_O.args))
        end
    end
    _O
end

function _simplify_constants(O)
    if Symbol(O.op) == :*
        # If any variable is `Constant(0)`, zero the whole thing
        # If any variable is `Constant(1)`, remove that `Constant(1)` unless
        # they are all `Constant(1)`, in which case simplify to a single variable
        if any(x->typeof(x)<:Variable && (isequal(x,Constant(0)) || isequal(x,Constant(-0))),O.args)
            return Constant(0)
        elseif any(x->typeof(x)<:Variable && isequal(x,Constant(1)),O.args)
            idxs = find(x->typeof(x)<:Variable && isequal(x,Constant(1)),O.args)
            _O = Operation(O.op,O.args[1:length(O.args) .∉ (idxs,)])
            if isempty(_O.args)
                return Constant(1)
            elseif length(_O.args) == 1
                return _O.args[1]
            else
                return _O
            end
        else
            return O
        end
    elseif Symbol(O.op) == :+ && any(x->typeof(x)<:Variable && (isequal(x,Constant(0)) || isequal(x,Constant(-0))),O.args)
        # If there are Constant(0)s in a big `+` expression, get rid of them
        idxs = find(x->typeof(x)<:Variable && (isequal(x,Constant(0)) || isequal(x,Constant(-0))),O.args)
        _O = Operation(O.op,O.args[1:length(O.args) .∉ (idxs,)])
        if isempty(_O.args)
            return Constant(0)
        elseif length(_O.args) == 1
            return _O.args[1]
        else
            return O
        end
    elseif O.op == identity
        return O.args[1]
    elseif Symbol(O.op) == :- && length(O.args) == 1
        return Operation(*,Expression[-1,O.args[1]])
    else
        return O
    end
end
simplify_constants(x::Variable) = x
_simplify_constants(x::Variable) = x

export simplify_constants

ext: jl


text:
abstract type RpTAlgorithm end

function get_rpt_alg()
    implemented = [ijk(), ijk2(), abc()]
    N = Options.get("pt_alg")
    try
        return implemented[N]
    catch BoundsError
        throw(FermiException("implementation number $N not available for RCCSD(T)."))
    end
end
"""
    Fermi.CoupledCluster.RCCSDpT

Fermi struct that holds information about RCCSD(T) wavefunctions

_struct tree:_

**RCCSD(T)** <: AbstractCCWavefunction <: AbstractCorrelatedWavefunction <: AbstractWavefunction
"""
struct RCCSDpT{T} <: AbstractCCWavefunction
    CCSD::RCCSD{T}
    energy::T
    correction::T
end

"""
    Fermi.CoupledCluster.get_pt_alg

Returns a singleton type corresponding to a RCCSD(T) implementation based on the options.
"""
function get_pt_alg(N::Int = Options.get("pt_alg"))
    try
        return get_pt_alg(Val(N))
    catch MethodError
        throw(FermiException("implementation number $N not available for RCCSD(T)."))
    end
end

# Implementations
struct ijk <: RpTAlgorithm end
include("ijk.jl")
get_pt_alg(x::Val{1}) = ijk()
struct ijk2 <: RpTAlgorithm end
include("ijk2.jl")
get_pt_alg(x::Val{2}) = ijk2()
struct abc <: RpTAlgorithm end
include("abc.jl")
get_pt_alg(x::Val{3}) = abc()

function RCCSDpT(x...)
    if !any(i-> i isa RpTAlgorithm, x)
        RCCSDpT(x..., get_rpt_alg())
    else
        # Print the type of arguments given for a better feedback
        args = "("
        for a in x[1:end-1]
            args *= "$(typeof(a)), "
        end
        args = args[1:end-2]*")"
        throw(FermiException("invalid arguments for RCCSD(T) method: $args"))
    end
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
abstract type AbstractDataEnvelopment end

# Abstract DEA model
abstract type AbstractDEA{S<:AbstractDataEnvelopment,T<:RS} <: AbstractArray{Any,1} end

function getdata(A::AbstractDEA)
	# Concrete types should overload this method to get array-like indexing features
end

# Return parametric types as a tuple
Base.eltype{S,T}(::Type{AbstractDEA{S,T}}) = (S, T)

Base.size(A::AbstractDEA) = size(getdata(A))
Base.IndexStyle(::Type{AbstractDEA}) = IndexLinear()

function Base.getindex(A::AbstractDEA, I)
	all(1 .<= I .<= size(A,1)) || throw(BoundsError(A, I))
	return getdata(A)[I]
end

# Broken since Julia 0.5! Will it ever come back?
# function Base.call{T<:AbstractDEA}(A::T)
# 	Data = getdata(A)
# 	res = Array(DEAResult,getnrdmu(Data))
#
# 	@sync @parallel for k in eachindex(Data)
# 		res[k] = Base.call(A,Data[k]...)
# 	end
#
# 	return res
# end

max_stars_repo_name: kepiej/NEPA


text:
Array{T}(::UndefInitializer, ::Tuple{PosInfinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{Integer, PosInfinity}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{PosInfinity, Integer}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{Integer, PosInfinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{PosInfinity, Integer}) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{PosInfinity, PosInfinity}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{PosInfinity, PosInfinity}) where T = throw(ArgumentError("Cannot create infinite Array"))

Array{T}(::UndefInitializer, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Array{T}(::UndefInitializer, ::PosInfinity, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::PosInfinity, ::Integer) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Integer, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Matrix{T}(::UndefInitializer, ::PosInfinity, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::PosInfinity, ::Integer) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Integer, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Vector{T}(::UndefInitializer, ::Tuple{PosInfinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Vector{T}(::UndefInitializer, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))

similar(A::AbstractArray, ::Type{T}, axes::Tuple{OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(A::AbstractArray, ::Type{T}, axes::Tuple{OneToInf{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{PosInfinity}) where T = cache(Zeros{T}(dims))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{PosInfinity,PosInfinity}) where T = cache(Zeros{T}(dims))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{Integer,PosInfinity}) where T = cache(Zeros{T}(dims))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{PosInfinity,Integer}) where T = cache(Zeros{T}(dims))

similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int},OneTo{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneTo{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{PosInfinity}) where T = cache(Zeros{T}(dims))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{PosInfinity,PosInfinity}) where T = cache(Zeros{T}(dims))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{Integer,PosInfinity}) where T = cache(Zeros{T}(dims))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{PosInfinity,Integer}) where T = cache(Zeros{T}(dims))

similar(arr::AbstractArray{T}, sz::Union{Infinity,Integer,AbstractUnitRange}...) where T = similar(arr, Base.to_shape(sz)...)
similar(arr::AbstractArray, ::Type{T}, sz::Union{Infinity,Integer,AbstractUnitRange}...) where T = similar(arr, T, Base.to_shape(sz)...)

for Typ in (:Zeros, :Ones)
    @eval begin
        @inline $Typ{T, N}(sz::Tuple{Vararg{Union{Infinity,Integer}, N}}) where {T, N} = $Typ{T,N}(oneto.(sz))
        @inline $Typ{T, N}(sz::Vararg{Union{Infinity,Integer}, N}) where {T, N} = $Typ{T,N}(sz)
        @inline $Typ{T}(sz::Vararg{Union{Infinity,Integer},N}) where {T, N} = $Typ{T, N}(sz)
        @inline $Typ{T}(sz::Infinity) where T = $Typ{T,1}(sz)
    end
end

Eye{T}(::Infinity) where T = Eye{T}(ℵ₀)
Eye{T}(::Infinity, ::Infinity) where T = Eye{T}(ℵ₀, ℵ₀)
Eye{T}(::Infinity, n::Integer) where T = Eye{T}(ℵ₀, n)
Eye{T}(m::Integer, ::Infinity) where T = Eye{T}(m, ℵ₀)
Eye(::Infinity) = Eye(ℵ₀)
Eye(::Infinity, ::Infinity) = Eye(ℵ₀, ℵ₀)
Eye(::Infinity, n::Integer) = Eye(ℵ₀, n)
Eye(m::Integer, ::Infinity) = Eye(m, ℵ₀)

@inline Fill(x, sz::Union{Infinity,Integer}...) = Fill(x, Base.to_shape(sz)...)
@inline Fill{T}(x, sz::Union{Infinity,Integer}...) where T = Fill{T}(x, Base.to_shape(sz)...)

zeros(sz::Union{Infinity, Integer, AbstractUnitRange}...) = zeros(Base.to_shape(sz)...)
zeros(::Type{T}, sz::Union{Infinity, Integer, AbstractUnitRange}...) where T = zeros(T, Base.to_shape(sz)...)
zeros(::Type{T}, ::Tuple{PosInfinity}) where T = cache(Zeros{T}(∞))
zeros(::Type{T}, nm::Tuple{Integer, PosInfinity}) where T = cache(Zeros{T}(nm...))
zeros(::Type{T}, nm::Tuple{PosInfinity, Integer}) where T = cache(Zeros{T}(nm...))
zeros(::Type{T}, nm::Tuple{PosInfinity, PosInfinity}) where T = cache(Zeros{T}(nm...))

ones(sz::Union{Infinity, Integer, AbstractUnitRange}...) = ones(Base.to_shape(sz)...)
ones(::Type{T}, ::Tuple{PosInfinity}) where T = cache(Ones{T}(∞))
ones(::Type{T}, nm::Tuple{Integer, PosInfinity}) where T = cache(Ones{T}(nm...))
ones(::Type{T}, nm::Tuple{PosInfinity, Integer}) where T = cache(Ones{T}(nm...))
ones(::Type{T}, nm::Tuple{PosInfinity, PosInfinity}) where T = cache(Ones{T}(nm...))

fill(x, sz::Union{Infinity, Integer, AbstractUnitRange}...) = fill(x, Base.to_shape(sz)...)
fill(x, ::Tuple{PosInfinity}) = cache(Fill(x,∞))
fill(x, nm::Tuple{Integer, PosInfinity}) = cache(Fill(x,nm...))
fill(x, nm::Tuple{PosInfinity, Integer}) = cache(Fill(x,nm...))
fill(x, nm::Tuple{PosInfinity, PosInfinity}) = cache(Fill(x,nm...))

# This gets called when infinit number of columns
axes_print_matrix_row(_, io, X, A, i, ::AbstractVector{<:PosInfinity}, sep) = nothing
print_matrix_row(io::IO, X::AbstractVecOrMat, A::Vector, i::Integer, cols::AbstractVector{<:PosInfinity}, sep::AbstractString) = nothing
Base.print_matrix_row(io::IO,
        X::Union{LayoutMatrix,
        LayoutVector,
        AbstractTriangular{<:Any,<:LayoutMatrix},
        AdjOrTrans{<:Any,<:LayoutMatrix},
        AdjOrTrans{<:Any,<:LayoutVector},
        HermOrSym{<:Any,<:LayoutMatrix},
        SubArray{<:Any,2,<:LayoutMatrix},
        Diagonal{<:Any,<:LayoutVector}}, A::Vector,
        i::Integer, cols::AbstractVector{<:PosInfinity}, sep::AbstractString) =
    axes_print_matrix_row(axes(X), io, X, A, i, cols, sep)
Base.print_matrix_row(io::IO,
        X::Union{AbstractFill{<:Any,1},
                 AbstractFill{<:Any,2},
                 Diagonal{<:Any,<:AbstractFill{<:Any,1}},
                 RectDiagonal,
                 AbstractTriangular{<:Any,<:AbstractFill{<:Any,2}}
                 }, A::Vector,
        i::Integer, cols::AbstractVector{<:PosInfinity}, sep::AbstractString) =
    axes_print_matrix_row(axes(X), io, X, A, i, cols, sep)

print_matrix_vdots(io::IO, vdots::AbstractString,
        A::Vector, sep::AbstractString, M::Integer, ::NotANumber, pad_right::Bool = true) = nothing

# Avoid infinite loops on maximum
Base.mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ::PosInfinity) =
    throw(ArgumentError("Cannot call mapreduce on an infinite length $(typeof(A))"))

function show_delim_array(io::IO, itr::AbstractArray, op, delim, cl,
                          delim_one, i1, ::PosInfinity)
    print(io, op)
    l = 20
    if !show_circular(io, itr)
        recur_io = IOContext(io, :SHOWN_SET => itr)
        if !haskey(io, :compact)
          recur_io = IOContext(recur_io, :compact => true)
        end
        first = true
        i = i1
        if 20 >= i1
          while true
              if !isassigned(itr, i)
                  print(io, undef_ref_str)
              else
                  x = itr[i]
                  show(recur_io, x)
              end
              i += 1
              if i > l
                  print(io, delim)
                  print(io, ' ')
                  print(io, '…')
                  delim_one && first && print(io, delim)
                  break
              end
              first = false
              print(io, delim)
              print(io, ' ')
          end
        end
    end
    print(io, cl)
end

#####
# FillArrays
#####

# Lazy Broadacasting
for typ in (:Ones, :Zeros, :Fill)
    @eval begin
        BroadcastStyle(::Type{$typ{T,N,NTuple{N,<:OneToInf}}}) where {T,N} = LazyArrayStyle{N}()
        BroadcastStyle(::Type{$typ{T,2,<:Tuple{<:Any,<:OneToInf}}}) where {T} = LazyArrayStyle{2}()
        BroadcastStyle(::Type{$typ{T,2,<:Tuple{<:OneToInf,<:Any}}}) where {T} = LazyArrayStyle{2}()
    end
end

BroadcastStyle(::Type{<:Diagonal{T,<:AbstractFill{T,1,Tuple{OneToInf{I}}}}}) where {T,I} = LazyArrayStyle{2}()

## Support broadcast(*, ::AbstractFill, A)

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Ones{T,1,Tuple{OneToInf{Int}}},AbstractArray{V,N}}}) where {N,T,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    convert(AbstractArray{promote_type(T,V),N}, b)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},Ones{V,1,Tuple{OneToInf{Int}}}}}) where {N,T,V}
    a,b = bc.args
    @assert bc.axes == axes(a)
    convert(AbstractArray{promote_type(T,V),N}, a)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractFill{T,1,Tuple{OneToInf{Int}}},AbstractArray{V,N}}}) where {N,T<:Number,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    getindex_value(a) * b
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},AbstractFill{V,1,Tuple{OneToInf{Int}}}}}) where {N,T,V<:Number}
    a,b = bc.args
    @assert bc.axes == axes(a)
    a * getindex_value(b)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractFill{T,1,Tuple{OneToInf{Int}}},AbstractArray{V,N}}}) where {N,T,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    Ref(getindex_value(a)) .* b # Use broadcast in-case a is array-valued
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},AbstractFill{V,1,Tuple{OneToInf{Int}}}}}) where {N,T,V}
    a,b = bc.args
    @assert bc.axes == axes(a)
    a .* Ref(getindex_value(b)) # Use broadcast in-case b is array-valued
end

# row Vector
function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Adjoint{T,Ones{T,1,Tuple{OneToInf{Int}}}},AbstractMatrix{V}}}) where {T,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    convert(AbstractMatrix{promote_type(T,V)}, b)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractMatrix{T},Adjoint{V,Ones{V,1,Tuple{OneToInf{Int}}}}}}) where {T,V}
    a,b = bc.args
    @assert bc.axes == axes(a)
    convert(AbstractMatrix{promote_type(T,V)}, a)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Adjoint{T,<:AbstractFill{T,1,Tuple{OneToInf{Int}}}},AbstractMatrix{V}}}) where {T,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    getindex_value(a') * b
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractMatrix{T},Adjoint{V,<:AbstractFill{V,1,Tuple{OneToInf{Int}}}}}}) where {T,V}
    a,b = bc.args
    @assert bc.axes == axes(a)
    a * getindex_value(b')
end

#####
# Diagonal
#####

BroadcastStyle(::Type{<:Diagonal{<:Any,<:AbstractInfUnitRange}}) = LazyArrayStyle{2}()

#####
# Vcat length
#####

function getindex(f::Vcat{T,1}, k::PosInfinity) where T
    length(f) == ℵ₀ || throw(BoundsError(f,k))
    ∞
end

_gettail(k, a::Number, b...) = k ≤ 1 ? tuple(a, b...) : _gettail(k - length(a), b...)
_gettail(k, a, b...) = k ≤ length(a) ? tuple(a[k:end], b...) : _gettail(k - length(a), b...)
_vcat(a) = a
_vcat(a, b, c...) = Vcat(a, b, c...)
getindex(A::Vcat, r::InfUnitRange) = Base.invoke(getindex, Tuple{AbstractArray, Any}, A, r)
_unsafe_getindex(::IndexLinear, A::Vcat, r::InfUnitRange) = _vcat(_gettail(first(r), A.args...)...)

getindex(A::AbstractVector, r::InfRanges{<:Integer}) = layout_getindex(A, r)
getindex(A::LayoutVector, r::InfRanges{<:Integer}) = layout_getindex(A, r)
getindex(A::CachedVector, r::InfRanges{<:Integer}) = layout_getindex(A, r)
getindex(A::AbstractFill{<:Any,1}, r::InfRanges{<:Integer}) = FillArrays.fillsimilar(A, length(r))

Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, I::AbstractFill) = Base.checkindex(Bool, inds, getindex_value(I))
LazyArrays.cache_getindex(::InfiniteCardinal{0}, A::AbstractVector, I, J...) = layout_getindex(A, I, J...)


text:
module TorchHub
using PyCallChainRules.Torch: TorchModuleWrapper, torch, functorch, ispysetup

using Statistics
using BenchmarkTools

using Zygote
using PyCallChainRules.Torch.PyCall

using CUDA

if !ispysetup[]
    return
end

py"""
import torch
def bn2group(module):
    num_groups = 16 # hyper_parameter of GroupNorm
    module_output = module
    if isinstance(module, torch.nn.modules.batchnorm._BatchNorm):
        module_output = torch.nn.GroupNorm(num_groups,
                                           module.num_features,
                                           module.eps,
                                           module.affine,
                                          )
        if module.affine:
            with torch.no_grad():
                module_output.weight = module.weight
                module_output.bias = module.bias
        module_output.running_mean = module.running_mean
        module_output.running_var = module.running_var
        module_output.num_batches_tracked = module.num_batches_tracked
        if hasattr(module, "qconfig"):
            module_output.qconfig = module.qconfig

    for name, child in module.named_children():
        module_output.add_module(name, bn2group(child))

    del module
    return module_output
"""

include("pytorch_utils.jl")

batchsizes = [1, ]

suite = BenchmarkGroup()

for device in ["cpu", "cuda"]
    ss = suite["$device"] = BenchmarkGroup()
    for batchsize in batchsizes
        model = torch.hub.load("pytorch/vision", "resnet18", pretrained=true)
        model_gn = py"bn2group"(model)

        s = ss["bs=$batchsize"] = BenchmarkGroup()
        # Forward pass
        s["forward"] = BenchmarkGroup()
        forward_pass!(s["forward"], model_gn, randn(Float32, 3, 224, 224), batchsize, device; samples=20)

        # Gradient pass
        s["backward"] = BenchmarkGroup()
        b2 = backward_pass!(s["backward"], model_gn, randn(Float32, 3, 224, 224), batchsize, device; samples=20)
    end
end
end
TorchHub.suite


text:
module Iteration

import Base.start
import Base.next
import Base.done
import Base.length

immutable Slicerator{T}
    array :: Array{T}
    length :: Uint64
    slice_size :: Uint64
end

function slicerator{T}(array :: Array{T}, slice_size)
    len = length(array)
    @assert(len % slice_size == 0)
    Slicerator(array, uint64(len), uint64(slice_size))
end

length(si :: Slicerator) = uint64(si.length / si.slice_size)
start(si :: Slicerator) = int64(0)
done(si :: Slicerator, i :: Int64) = (i * si.slice_size) >= length(si.array)
function next(si :: Slicerator, i :: Int64)
    sidx = int64((i * si.slice_size) + 1)
    eidx = int64((i+1) * si.slice_size)

    (slice(si.array, sidx:eidx), i+1)
end

immutable OffsetIterator
    inner
    offset :: Int64
end

offset(offset :: Int64, inner) = OffsetIterator(inner, offset)

function length(oi :: OffsetIterator)
    len = length(oi.inner)
    if len > oi.offset
        len - oi.offset
    else
        0
    end
end
start(oi :: OffsetIterator) = start(oi.inner) + oi.offset
done(oi :: OffsetIterator, i) = done(oi.inner, i)
next(oi :: OffsetIterator, i) = next(oi.inner, i)

export slicerator, offset

end # module Iteration


text:
module FluxMO

    # package code goes here
    include("Clustering.jl")
    include("Validation.jl")

end # module

max_stars_repo_name: swiesend/FluxMO

max_issues_repo_name: swiesend/FluxMO

max_forks_repo_name: swiesend/FluxMO


text:
#=
This file contains functional modules (modules without any learnable contents).
Functional Modules don't take ctx as an argument
=#

abstract type FnModule <: KnetModule end

# Activation
abstract type Activation <: FnModule end

for (Act, fn) in zip([:ReLU, :Sigmoid, :Tanh], [:relu, :sigm, :tanh])
    eval(:(
        begin
           type $(Act) <: Activation; ;end
           (m::$(Act))(x) = $(fn).(x)
           #forward(ctx, m::$(Act), x) = $(fn).(x)
        end
    ))
end

abstract type Pool <: FnModule end

# Pooling
type MaxPool <: Pool; opt; end

MaxPool(;o...) = MaxPool(o)

(m::MaxPool)(x) = pool(x; m.opt..., mode=0)

type AvgPool <: Pool
    opt
    include_pad::Bool
end

AvgPool(;include_pad=true, o...) = AvgPool(o, include_pad)

(m::AvgPool)(x) =
    pool(x; m.opt..., mode=1+Int(m.include_pad))

# Dropout
type Dropout <: FnModule
    p::AbstractFloat
    train::Union{Bool, Void}
end

Dropout(pdrop; train=true) = Dropout(pdrop, train)

function (d::Dropout)(x)
    train = d.train
    if ~isa(train, Bool)
        train = isa(x, AutoGrad.Rec)
    end
    p = train ? d.p : 0
    return dropout(x, p)
end


text:
#=
run native instance tests from test/nativeinstances.jl and display basic benchmarks
=#

using Test
using Printf
import DataFrames
import Hypatia
import Hypatia.Solvers
include(joinpath(@__DIR__, "nativeinstances.jl"))
include(joinpath(@__DIR__, "nativesets.jl"))

# default options to solvers
default_options = (
    # verbose = true,
    verbose = false,
    default_tol_relax = 10,
    )

all_reals = [
    Float32,
    Float64,
    BigFloat,
    ]
diff_reals = [
    Float64,
    BigFloat,
    ]

string_nameof(T) = string(nameof(T))

type_name(::T) where T = string_nameof(T)

function test_instance_solver(
    inst_name::String,
    T::Type{<:Real},
    options::NamedTuple,
    test_info::String = "",
    )
    test_info = "$inst_name $T $test_info"
    @testset "$test_info" begin
        println(test_info, " ...")
        solver = Solvers.Solver{T}(; options...)
        test_time = @elapsed eval(Symbol(inst_name))(T, solver = solver)
        push!(perf, (inst_name, string(T), type_name(solver.stepper),
            type_name(solver.syssolver), solver.init_use_indirect,
            solver.preprocess, solver.reduce, test_time,
            string(Solvers.get_status(solver))))
        @printf("%8.2e seconds\n", test_time)
    end
    return nothing
end

perf = DataFrames.DataFrame(
    inst_name = String[],
    real_T = String[],
    stepper = String[],
    syssolver = String[],
    init_use_indirect = Bool[],
    preprocess = Bool[],
    reduce = Bool[],
    test_time = Float64[],
    status = String[],
    )

@testset "native tests" begin

@testset "default options tests" begin
    println("starting default options tests")
    inst_defaults = vcat(
        inst_preproc,
        inst_infeas,
        inst_cones_many,
        )
    for inst_name in inst_defaults
        test_instance_solver(inst_name, Float64, default_options)
    end
end

@testset "no preprocess tests" begin
    println("\nstarting no preprocess tests")
    for inst_name in inst_cones_few, T in diff_reals
        options = (; default_options..., preprocess = false, reduce = false,
            syssolver = Solvers.SymIndefDenseSystemSolver{T}())
        test_instance_solver(inst_name, T, options)
    end
end

@testset "indirect solvers tests" begin
    println("\nstarting indirect solvers tests")
    for inst_name in inst_indirect, T in diff_reals
        options = (; default_options..., init_use_indirect = true,
            preprocess = false, reduce = false,
            syssolver = Solvers.SymIndefIndirectSystemSolver{T}(),
            tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4,
            tol_infeas = 1e-6)
        test_instance_solver(inst_name, T, options)
    end
end

@testset "system solvers tests" begin
    println("\nstarting system solvers tests")
    syssolvers = [
        (Solvers.NaiveDenseSystemSolver, diff_reals),
        (Solvers.NaiveSparseSystemSolver, [Float64,]),
        (Solvers.NaiveElimDenseSystemSolver, diff_reals),
        (Solvers.NaiveElimSparseSystemSolver, [Float64,]),
        (Solvers.SymIndefDenseSystemSolver, all_reals),
        (Solvers.SymIndefSparseSystemSolver, [Float64,]),
        (Solvers.QRCholDenseSystemSolver, all_reals),
        ]
    for inst_name in inst_minimal, (syssolver, real_types) in syssolvers,
        T in real_types
        options = (; default_options..., syssolver = syssolver{T}(),
            reduce = false)
        test_instance_solver(inst_name, T, options, string_nameof(syssolver))
    end
end

@testset "PredOrCentStepper tests" begin
    verbose = true
    println("\nstarting PredOrCentStepper tests (with printing)")

    # adjustment and curve search
    use_adj_curv = [(false, false), (true, false), (true, true)]
    for inst_name in inst_minimal, (adj, curv) in use_adj_curv, T in diff_reals
        stepper = Solvers.PredOrCentStepper{T}(;
            use_adjustment = adj, use_curve_search = curv)
        options = (; default_options..., verbose = verbose, stepper = stepper)
        test_instance_solver(inst_name, T, options, "adj=$adj curv=$curv")
    end

    # other options
    for inst_name in inst_minimal
        T = Float64
        stepper = Solvers.PredOrCentStepper{T}(;
            # stepper options
            use_adjustment = false, use_curve_search = false,
            max_cent_steps = 8, pred_prox_bound = 0.0332,
            # searcher options
            min_prox = 0.0, prox_bound = 0.2844, use_max_prox = false,
            alpha_sched = [0.9999 * 0.7^i for i in 0:22])
        options = (; default_options..., verbose = verbose, stepper = stepper)
        test_instance_solver(inst_name, T, options, "other")
    end
end

@testset "CombinedStepper tests" begin
    verbose = true
    println("\nstarting CombinedStepper tests (with printing)")
    shifts = [0, 2]
    for inst_name in inst_minimal, shift in shifts, T in diff_reals
        options = (; default_options..., verbose = verbose,
            stepper = Solvers.CombinedStepper{T}(shift_sched = shift))
        test_instance_solver(inst_name, T, options, "shift=$shift")
    end
end

# println("\n")
# DataFrames.show(perf, allrows = true, allcols = true)
# println("\n")
end
;


text:
export FiniteException, write_solutions, read_solutions, write_parameters, read_parameters

"""
    write_solutions(filename, solutions)

Stores the as a plain text file onto disk.
The storage format is as follows. The first line indicates the number of solutions stored
followed by a blank line.
Then the solutions are stored where each solution is separated by a blank line.
Note that the solutions are *always* considered as complex numbers.
See [`read_solutions`](@ref) for reading the solution file back.

Note that this is the same file format as used in [Bertini](https://bertini.nd.edu).

## Example

```julia
julia> write_solutions("solutions.txt", [[1, 1], [-1, 2]])

shell> cat solutions.txt
2

1 0
1 0

-1 0
2 0

julia> read_solutions("solutions.txt")
2-element Array{Array{Complex{Int64},1},1}:
 [1 + 0im, 1 + 0im]
 [-1 + 0im, 2 + 0im]
```
"""
function write_solutions(filename, S::AbstractVector{<:AbstractVector{<:Number}})
    open(filename, "w") do f
        write(f, "$(length(S))\n")
        for s in S
            write(f, "\n")
            DelimitedFiles.writedlm(f, [real.(s) imag.(s)], ' ')
        end
    end
end

"""
    read_solutions(filename)

Read the solutions stored with [`write_solutions`](@ref).

## Example

```julia
julia> write_solutions("solutions.txt", [[1, 1], [-1, 2]])

julia> read_solutions("solutions.txt")
2-element Array{Array{Complex{Int64},1},1}:
 [1 + 0im, 1 + 0im]
 [-1 + 0im, 2 + 0im]
```
"""
function read_solutions(filename)
    A = DelimitedFiles.readdlm(filename; skipblanks = false)
    n = convert(Int, A[1, 1])
    # figure out number of variables
    nvars = 0
    for k = 4:size(A, 1)
        if isempty(A[k, 1])
            nvars = k - 3
            break
        end
    end
    map(3:nvars+1:size(A, 1)) do i
        map(i:i+nvars-1) do k
            A[k, 1] + im * A[k, 2]

        end
    end
end

"""
    write_parameters(filename, parameters)

Stores the parameters as a plain text file onto disk.
The storage format is as follows. The first line indicates the number of parameter values
stored followed by a blank line. Then the parameter values are stored.

See [`read_parameters`](@ref)

Note that this is the same file format as used in [Bertini](https://bertini.nd.edu).

## Example

```julia
julia> write_parameters("parameters.txt", [2.0, -3.2 + 2im])

shell> cat parameters.txt
2

2.0 0.0
-3.2 2.0

julia> read_parameters("parameters.txt")
2-element Array{Complex{Float64},1}:
  2.0 + 0.0im
 -3.2 + 2.0im
```
"""
function write_parameters(filename, p::Vector{<:Number})
    open(filename, "w") do f
        write(f, "$(length(p))\n\n")
        DelimitedFiles.writedlm(f, [real.(p) imag.(p)], ' ')
    end
end

"""
    read_parameters(filename)

Read the parameters stored with [`write_parameters`](@ref).

## Example

```julia
julia> write_parameters("parameters.txt", [2.0, -3.2 + 2im])

julia> read_parameters("parameters.txt")
2-element Array{Complex{Float64},1}:
  2.0 + 0.0im
 -3.2 + 2.0im
```
"""
function read_parameters(filename)
    A = DelimitedFiles.readdlm(filename)
    n = A[1, 1]
    map(1:n) do i
        A[i+1, 1] + im * A[i+1, 2]
    end
end

###
### Exceptions
###
struct KeywordArgumentException <: Exception
    key::Symbol
    given::Any
    msg::String
end
function KeywordArgumentException(key, given)
    KeywordArgumentException(key, given, "")
end
function Base.showerror(io::IO, E::KeywordArgumentException)
    print(io, "Invalid argument $given for $key. ", msg)
end

struct FiniteException <: Exception
    dim::Int
end
function Base.showerror(io::IO, E::FiniteException)
    print(
        io,
        "FiniteException: The solution set of the given system has at least dimension " *
        "$(E.dim) > 0. Consider intersecting with an (affine) subspace of codimension " *
        "$(E.dim) to reduce to (possibly) finitely many solutions.",
    )
end

@noinline function unsupported_kwargs(kwargs)
    if !(isempty(kwargs))
        msg = join(["$k = $v" for (k, v) in pairs(kwargs)], ", ")
        @warn "Ingored unsupported keyword arguments: $msg"
    end
end

# Perform type piracy here until https://github.com/JuliaArrays/StructArrays.jl/issues/131
# is fixed
const _C64SA = StructArrays.StructArray{
    Complex{Float64},
    2,
    NamedTuple{(:re, :im),Tuple{Array{Float64,2},Array{Float64,2}}},
    Int64,
}
Base.@propagate_inbounds function Base.setindex!(s::_C64SA, vals, I::Int)
    @boundscheck checkbounds(s, I)
    StructArrays.foreachfield(
        (col, val) -> (@inbounds col[I] = val),
        s,
        convert(ComplexF64, vals),
    )
    s
end

"""
    rand_unitary_matrix(n::Int, T=ComplexF64)

Samples a `n × n` unitary Matrix uniformly from the space of all unitary n × n matrices.

See https://arxiv.org/abs/math-ph/0609050 for a derivation.
"""
function rand_unitary_matrix(n::Int, T::Type = ComplexF64)
    Z = randn(T, n, n) ./ sqrt(2)
    Q, R = LA.qr(Z)
    Λ = LA.diagm(0 => [R[i, i] / abs(R[i, i]) for i = 1:n])
    Q * Λ
end

fast_abs(z::Complex) = sqrt(abs2(z))
fast_abs(x::Real) = abs(x)

"Like `min(a,b)`` but ignoring any `NaN` values."
nanmin(a, b) = isnan(a) ? b : (isnan(b) ? a : min(a, b))
nanmax(a, b) = isnan(a) ? b : (isnan(b) ? a : max(a, b))

always_false(::Any) = false

"""
    all2(f, a, b)

Checks if `f(a[i], b[i]) == true` for `i=1:length(a)`.
"""
function all2(f::F, a::AbstractVector, b::AbstractVector) where {F}
    length(a) == length(b) || return false
    for (ai, bi) in zip(a, b)
        f(ai, bi) || return false
    end
    true
end

"""
    unpack(a::Union{Nothing, T}, b::T)

Returns `a` if it is not `nothing`, otherwise `b`.
"""
unpack(a::Union{Nothing,T}, b::T) where {T} = a === nothing ? b : a

plural(singularstr, n) = n == 1 ? singularstr : singularstr * "s"

"""
     print_fieldnames(io::IO, obj)

 A better default printing for structs.
 """
function print_fieldnames(io::IO, obj)
    println(io, typeof(obj), ":")
    for name in fieldnames(typeof(obj))
        print_fieldname(io, obj, name)
    end
end
function print_fieldname(io::IO, obj, name)
    if getfield(obj, name) !== nothing
        val = getfield(obj, name)
        print(io, " • ", name, " → ")
        if val isa AbstractFloat
            println(io, round(val; sigdigits = 5))
        else
            println(io, val)
        end
    end
end

"""
    @tspawnat tid -> task
Mimics `Base.Threads.@spawn`, but assigns the task to thread `tid`.
# Example
```julia
julia> t = @tspawnat 4 Threads.threadid()
Task (runnable) @0x0000000010743c70
julia> fetch(t)
4
```
"""
macro tspawnat(thrdid, expr)
    thunk = esc(:(() -> ($expr)))
    var = esc(Base.sync_varname)
    tid = esc(thrdid)
    quote
        if $tid < 1 || $tid > Threads.nthreads()
            throw(
                AssertionError(
                    "@tspawnat thread assignment ($($tid)) must be between 1 and Threads.nthreads() (1:$(Threads.nthreads()))",
                ),
            )
        end
        local task = Task($thunk)
        task.sticky = false
        ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, $tid - 1)
        if $(Expr(:isdefined, var))
            push!($var, task)
        end
        schedule(task)
        task
    end
end

mutable struct SegmentStepper
    start::ComplexF64
    target::ComplexF64
    abs_Δ::Float64
    forward::Bool
    # current
    s::Float64
    # proposed
    s′::Float64
end
SegmentStepper(start::Number, target::Number) =
    SegmentStepper(ComplexF64(start), ComplexF64(target))
SegmentStepper(start::ComplexF64, target::ComplexF64) =
    init!(SegmentStepper(start, target, 0.0, true, 0.0, 0.0), start, target)

init!(S::SegmentStepper, start::Number, target::Number) =
    init!(S, ComplexF64(start), ComplexF64(target))
function init!(S::SegmentStepper, start::ComplexF64, target::ComplexF64)
    S.start = start
    S.target = target
    S.abs_Δ = abs(target - start)
    S.forward = abs(start) < abs(target)
    S.s = S.s′ = S.forward ? 0.0 : S.abs_Δ
    S
end

is_done(S::SegmentStepper) = S.forward ? S.s == S.abs_Δ : S.s == 0.0

step_success!(S::SegmentStepper) = (S.s = S.s′; S)
function propose_step!(S::SegmentStepper, Δs::Real)
    if S.forward
        S.s′ = min(S.s + Δs, S.abs_Δ)
    else
        S.s′ = max(S.s - Δs, 0.0)
    end
    S
end
dist_to_target(S::SegmentStepper) = S.forward ? S.abs_Δ - S.s : S.s

function Base.getproperty(S::SegmentStepper, sym::Symbol)
    if sym == :Δs
        s = getfield(S, :s)
        s′ = getfield(S, :s′)
        Δs = getfield(S, :forward) ? s′ - s : s - s′
        return Δs
    elseif sym === :t
        s = getfield(S, :s)
        start = getfield(S, :start)
        target = getfield(S, :target)
        abs_Δ = getfield(S, :abs_Δ)
        forward = getfield(S, :forward)
        return _t_helper(start, target, s, abs_Δ, forward)
    elseif sym == :t′
        s′ = getfield(S, :s′)
        start = getfield(S, :start)
        target = getfield(S, :target)
        abs_Δ = getfield(S, :abs_Δ)
        forward = getfield(S, :forward)
        return _t_helper(start, target, s′, abs_Δ, forward)
    elseif sym == :Δt
        s = getfield(S, :s)
        s′ = getfield(S, :s′)
        start = getfield(S, :start)
        target = getfield(S, :target)
        abs_Δ = getfield(S, :abs_Δ)
        forward = getfield(S, :forward)
        if forward
            Δt = ((s′ - s) / abs_Δ) * (target - start)
        else
            Δt = ((s - s′) / abs_Δ) * (target - start)
        end
        return Δt
    else # fallback to getfield
        return getfield(S, sym)
    end
end
function _t_helper(start, target, s, Δ, forward)
    if forward
        if s == 0.0
            return start
        elseif s == Δ
            return target
        else
            return start + (s / Δ) * (target - start)
        end
    else
        if s == Δ
            return start
        elseif s == 0.0
            return target
        else
            return target + (s / Δ) * (start - target)
        end
    end
end

function Base.show(io::IO, val::SegmentStepper)
    print(io, "SegmentStepper:")
    for field in [:start, :target, :t, :Δt]
        print(io, "\n • ", field, " → ", getproperty(val, field))
    end
end

"""
    nthroot(x, n)

Compute the `n`-th root of `x`.
"""
function nthroot(x::Real, N::Integer)
    if N == 4
        √(√(x))
    elseif N == 2
        √(x)
    elseif N == 3
        cbrt(x)
    elseif N == 1
        x
    elseif N == 0
        one(x)
    else
        x^(1 / N)
    end
end

import Base: MPFR
import Base.MPFR: MPFRRoundingMode

function set!(x::BigFloat, y::BigFloat, r::MPFRRoundingMode = MPFR.ROUNDING_MODE[])
    ccall(
        (:mpfr_set, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode),
        x,
        y,
        r,
    )
    x
end

function set!(x::BigFloat, d::Float64, r::MPFRRoundingMode = MPFR.ROUNDING_MODE[])
    ccall(
        (:mpfr_set_d, :libmpfr),
        Int32,
        (Ref{BigFloat}, Float64, MPFRRoundingMode),
        x,
        d,
        r,
    )
    x
end

# Basic arithmetic without promotion
for (fJ, fC) in ((:add!, :add), (:mul!, :mul))
    @eval begin
        # BigFloat
        function ($fJ)(z::BigFloat, x::BigFloat, y::BigFloat)
            ccall(
                ($(string(:mpfr_, fC)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode),
                z,
                x,
                y,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end

        function ($fJ)(z::BigFloat, x::BigFloat, c::Float64)
            ccall(
                ($(string(:mpfr_, fC, :_d)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Cdouble, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        ($fJ)(c::Float64, x::BigFloat) = ($fJ)(x, c)

        function ($fJ)(z::BigFloat, x::BigFloat, c::Int64)
            ccall(
                ($(string(:mpfr_, fC, :_si)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Clong, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        ($fJ)(c::Int64, x::BigFloat) = ($fJ)(x, c)

        # BigInt
        function ($fJ)(z::BigFloat, x::BigFloat, c::BigInt)
            ccall(
                ($(string(:mpfr_, fC, :_z)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Ref{BigInt}, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        ($fJ)(c::BigInt, x::BigFloat) = ($fJ)(x, c)
    end
end

for (fJ, fC) in ((:sub!, :sub), (:div!, :div))
    @eval begin
        # BigFloat
        function ($fJ)(z::BigFloat, x::BigFloat, y::BigFloat)
            ccall(
                ($(string(:mpfr_, fC)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode),
                z,
                x,
                y,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end

        # Float32/Float64
        function ($fJ)(z::BigFloat, x::BigFloat, c::Float64)
            ccall(
                ($(string(:mpfr_, fC, :_d)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Cdouble, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        function ($fJ)(z::BigFloat, c::Float64, x::BigFloat)
            ccall(
                ($(string(:mpfr_, :d_, fC)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Cdouble, Ref{BigFloat}, MPFRRoundingMode),
                z,
                c,
                x,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end

        # BigInt
        function ($fJ)(z::BigFloat, x::BigFloat, c::BigInt)
            ccall(
                ($(string(:mpfr_, fC, :_z)), :libmpfr),
                Int32,
                (Ref{BigFloat}, Ref{BigFloat}, Ref{BigInt}, MPFRRoundingMode),
                z,
                x,
                c,
                MPFR.ROUNDING_MODE[],
            )
            return z
        end
        # no :mpfr_z_div function
    end
end

function sub!(z::BigFloat, c::BigInt, x::BigFloat)
    ccall(
        (:mpfr_z_sub, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigInt}, Ref{BigFloat}, MPFRRoundingMode),
        z,
        c,
        x,
        MPFR.ROUNDING_MODE[],
    )
    return z
end

function rem!(z::BigFloat, x::BigFloat, y::BigFloat, ::RoundingMode{:Nearest})
    ccall(
        (:mpfr_remainder, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode),
        z,
        x,
        y,
        MPFR.MPFRRoundNearest,
    )
    return z
end


text:
macro do_while(condition, block)
    quote
        let
            $block
            while $condition
                $block
            end
        end
    end |> esc
end

function _reg(s, quoted_attr, attr_name)
    @eval begin
        #$get_attr_name(x)      = getattr( x, $quoted_attr)
        #$set_attr_name(x, val) = setattr!(x, $quoted_attr, val)
        #export $get_attr_name
        #export $set_attr_name
        $attr_name(o::Object) = get_attr(o, $quoted_attr)
        $attr_name(w::World, id::Symbol) = get_attr(w, id, $quoted_attr)
        $attr_name(o::Object, v::Any) = set_attr!(o, $quoted_attr, v)
        $attr_name(w::World, id::Symbol, v::Any) = set_attr!(w, id, $quoted_attr, v)
        $attr_name(w::World, o::Object) = $attr_name(o)
        $attr_name(w::World, o::Object, v::Any) = $attr_name(o, v)
        export $attr_name
    end
end

macro register_attribute(attr)
    #get_attr_name = Symbol("get_", attr)
    #set_attr_name = Symbol("set_", attr)
    attr_name = Symbol(attr)
    quoted_attr  = QuoteNode(attr)
    s = string(attr)
    _reg(s, quoted_attr, attr_name)
end


text:
using NLSolvers, Test, DoubleFloats
function myfun(∇f, x)
    myfun(nothing, ∇f, x)
end
function myfun(∇²f, ∇f, x::T) where T
    if !(∇²f == nothing)
        ∇²f = 12x^2 - sin(x)
    end
    if !(∇f == nothing)
        ∇f = 4x^3 + cos(x)
    end

    fx = x^4 +sin(x)
    if ∇f == nothing && ∇²f == nothing
        return T(fx)
    elseif ∇²f == nothing
        return T(fx), T(∇f)
    else
        return T(fx), T(∇f), T(∇²f)
    end
end
@testset "scalar return types" begin
for T = (Float16, Float32, Float64, Rational{BigInt}, Double32, Double64)
    for M in (SR1, BFGS, DFP, Newton)
        if M == Newton
            obj = TwiceDiffed(myfun)
            res = minimize(obj, T(4), M())
            @test all(isa.(res[1:3], T))
        else
            obj = OnceDiffed(myfun)
            res = minimize(obj, T(4), M(Direct()))
            @test all(isa.(res[1:3], T))
            res = minimize(obj, T(4), M(Inverse()))
            @test all(isa.(res[1:3], T))
        end
    end
end
end


text:
#=
hypograph of geometric mean
(u, w) : wᵢ ≥ 0, u ≤ geom(w)
where geom(w) = (∏ᵢ wᵢ)^(1/d)

dual cone is
(u, w) : u ≤ 0, wᵢ ≥ 0, u ≥ -d * geom(w)

extended formulation
NOTE slightly different to EF in Hypatia solver paper
∃ θ ≥ u, ∃ λ, Σᵢ λᵢ ≥ 0, (λᵢ, θ, wᵢ) ∈ HypoPerLog
i.e. θ ≥ 0, wᵢ ≥ 0, λᵢ ≤ θ log(wᵢ / θ)
dual of HypoPerLog is (p, q, r) : p ≤ 0, r ≥ 0, q ≥ p * (log(-r / p) + 1)
=#

mutable struct HypoGeoMean{E <: NatExt} <: Cache
    oa_s::Vector{AE}
    d::Int
    θ::VR
    λ::Vector{VR}
    HypoGeoMean{E}() where {E <: NatExt} = new{E}()
end

function MOIPajarito.Cones.create_cache(
    oa_s::Vector{AE},
    cone::Hypatia.HypoGeoMeanCone{RealF},
    opt::Optimizer,
)
    @assert !cone.use_dual # TODO
    dim = MOI.dimension(cone)
    @assert dim == length(oa_s)
    d = dim - 1
    E = nat_or_ext(opt, d)
    cache = HypoGeoMean{E}()
    cache.oa_s = oa_s
    cache.d = d
    return cache
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::HypoGeoMean,
    opt::Optimizer,
)
    return _get_cuts(z[2:end], cache, opt)
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::HypoGeoMean,
    opt::Optimizer,
)
    us = s[1]
    @views ws = s[2:end]
    if us < opt.tol_feas || geomean(ws) - us > -opt.tol_feas
        return AE[]
    end

    # gradient cut is (-1, geom(w) / d ./ w)
    w_pos = max.(ws, 1e-7)
    c1 = geomean(w_pos) / cache.d
    r = c1 ./ w_pos
    return _get_cuts(r, cache, opt)
end

# unextended formulation

function MOIPajarito.Cones.add_init_cuts(cache::HypoGeoMean{Nat}, opt::Optimizer)
    # add variable bounds wᵢ ≥ 0
    @views w = cache.oa_s[2:end]
    JuMP.@constraint(opt.oa_model, w .>= 0)
    opt.use_init_fixed_oa || return

    # add cut (-d, e) on (u, w), a linearization at w = e
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, -cache.d * u + sum(w) >= 0)
    return
end

function _get_cuts(r::Vector{RealF}, cache::HypoGeoMean{Nat}, opt::Optimizer)
    # strengthened cut is (-d * geom(r), r)
    clean_array!(r) && return AE[]
    p = -cache.d * geomean(r)
    u = cache.oa_s[1]
    @views w = cache.oa_s[2:end]
    cut = JuMP.@expression(opt.oa_model, p * u + JuMP.dot(r, w))
    return [cut]
end

# extended formulation

MOIPajarito.Cones.num_ext_variables(cache::HypoGeoMean{Ext}) = 1 + cache.d

function MOIPajarito.Cones.extend_start(
    cache::HypoGeoMean{Ext},
    s_start::Vector{RealF},
    opt::Optimizer,
)
    @views w_start = s_start[2:end]
    if any(<(1e-9), w_start)
        return zeros(1 + cache.d)
    end
    θ_start = geomean(w_start)
    λ_start = [θ_start * log(w_i / θ_start) for w_i in w_start]
    return vcat(θ_start, λ_start)
end

function MOIPajarito.Cones.setup_auxiliary(cache::HypoGeoMean{Ext}, opt::Optimizer)
    @assert cache.d >= 2
    θ = cache.θ = JuMP.@variable(opt.oa_model, lower_bound = 0)
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, θ >= u)
    λ = cache.λ = JuMP.@variable(opt.oa_model, [1:(cache.d)])
    JuMP.@constraint(opt.oa_model, sum(λ) >= 0)
    return vcat(θ, λ)
end

function MOIPajarito.Cones.add_init_cuts(cache::HypoGeoMean{Ext}, opt::Optimizer)
    # add variable bounds wᵢ ≥ 0
    @views w = cache.oa_s[2:end]
    JuMP.@constraint(opt.oa_model, w .>= 0)
    opt.use_init_fixed_oa || return

    # add disaggregated cuts (-1, -1, 1) on (λᵢ, θ, wᵢ) TODO check
    θ = cache.θ
    λ = cache.λ
    JuMP.@constraint(opt.oa_model, [i in 1:(cache.d)], -λ[i] - θ + w[i] >= 0)
    return
end

function _get_cuts(r::Vector{RealF}, cache::HypoGeoMean{Ext}, opt::Optimizer)
    clean_array!(r) && return AE[]
    p = -geomean(r)
    @views w = cache.oa_s[2:end]
    θ = cache.θ
    λ = cache.λ
    cuts = AE[]
    for i in 1:(cache.d)
        r_i = r[i]
        iszero(r_i) && continue # TODO ?
        # strengthened disaggregated cut on (λᵢ, θ, wᵢ) is
        # p = -geom(r), (p, p * (log(-rᵢ / p) + 1), rᵢ)
        # TODO check math
        q = p * (log(-r_i / p) + 1)
        cut = JuMP.@expression(opt.oa_model, p * λ[i] + q * θ + r_i * w[i])
        push!(cuts, cut)
    end
    return cuts
end


text:
# # Plotting
#
# ## Basics
# There are a few different ways to plot in Julia (including calling PyPlot).

#
# Here we'll show you how to use `Plots.jl`.  If it's not installed yet, you need to use the package manager to install it, and Julia will precompile it for you the first time you use it:

## using Pkg
## Pkg.add("Plots")
using Plots

# One of the advantages to `Plots.jl` is that it allows you to seamlessly change backends. In this notebook, we'll try out the `gr()` and `plotlyjs()` backends.

#
# In the name of scientific inquiry, let's use this notebook to examine the relationship between the global temperature and the number of pirates between roughly 1860 and 2000.

globaltemperatures = [14.4, 14.5, 14.8, 15.2, 15.5, 15.8]
numpirates = [45000, 20000, 15000, 5000, 400, 17];

# Plots supports multiple backends — that is, libraries that actually do the drawing — all with the same API. To start out, let's try the GR backend.  You choose it with a call to `gr()`:

gr()

# and now we can use commands like `plot` and `scatter` to generate plots.

plot(numpirates, globaltemperatures, label="line")
scatter!(numpirates, globaltemperatures, label="points")

# The `!` at the end of the `scatter!` function name makes `scatter!` a mutating function, indicating that the scattered points will be added onto the pre-existing plot.
#
# In contrast, see what happens when you replace `scatter!` in the above with the non-mutating function `scatter`.
#
# Next, let's update this plot with the `xlabel!`, `ylabel!`, and `title!` commands to add more information to our plot.

xlabel!("Number of Pirates [Approximate]")
ylabel!("Global Temperature (C)")
title!("Influence of pirate population on global warming")

# This still doesn't look quite right. The number of pirates has decreased since 1860, so reading the plot from left to right is like looking backwards in time rather than forwards. Let's flip the x axis to better see how pirate populations have caused global temperatures to change over time!

xflip!()

# And there we have it!
#
# Note: We've had some confusion about this exercise. :) This is a joke about how people often conflate correlation and causation.
#
# **Without changing syntax, we can create this plot with the UnicodePlots backend**

Pkg.add("UnicodePlots")
unicodeplots()

#-

plot(numpirates, globaltemperatures, label="line")
scatter!(numpirates, globaltemperatures, label="points")
xlabel!("Number of Pirates [Approximate]")
ylabel!("Global Temperature (C)")
title!("Influence of pirate population on global warming")

# And notice how this second plot differs from the first!  Using text like this is a little silly in a Jupyter notebook where we have fancy drawing capabilities, but it can be very useful for quick and dirty visualization in a terminal.

#-

# ### Exercises
#
# #### 8.1
# Given
# ```julia
# x = -10:10
# ```
# plot y vs. x for $y = x^2$.  You may want to change backends back again.

# #### 8.2
# Execute the following code

p1 = plot(x, x)
p2 = plot(x, x.^2)
p3 = plot(x, x.^3)
p4 = plot(x, x.^4)
plot(p1, p2, p3, p4, layout = (2, 2), legend = false)

# and then create a $4x1$ plot that uses `p1`, `p2`, `p3`, and `p4` as subplots.

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
## Vector/Collections API

function append!(x::RLEVector, y::RLEVector)
    last_x_run = nrun(x)
    length_x = length(x)
    if last(x) == first(y)
        last_x_run = last_x_run - 1
    end
    new_nrun = last_x_run + nrun(y)
    resize!(x.runvalues, new_nrun)
    resize!(x.runends, new_nrun)
    x.runvalues[(last_x_run .+ 1:end)] = y.runvalues
    x.runends[(last_x_run .+ 1:end)] = y.runends .+ length_x
    x
end

function vcat(x::RLEVector, y::RLEVector...)
    out = copy(x)
    for yi in y
        append!(out, yi)
    end
    out
end

function pop!(x::RLEVector)
  runcount = nrun(x)
  isempty(x) && throw(ArgumentError("array must be non-empty"))
  item = x.runvalues[end]
  x.runends[end] -= 1
  if x.runends[end] == 0
    deleteat!(x.runvalues,runcount)
    deleteat!(x.runends,runcount)
  end
  return(item)
end

function push!(x::RLEVector{T1,T2},item) where {T1,T2 <: Integer}
  item = convert(T1,item) # Copying how base does it for arrays
  if !isempty(x) && (item == x.runvalues[end])
    x.runends[end] += 1
  else
    push!(x.runvalues,item)
    push!(x.runends, (isempty(x) ? 0 : x.runends[end]) + 1)
  end
  return(x)
end

function popfirst!(x::RLEVector)
  isempty(x) && throw(ArgumentError("array must be non-empty"))
  item = x.runvalues[1]
  x.runends[:] = x.runends .- 1
  if x.runends[1] == 0
    deleteat!(x.runvalues,1)
    deleteat!(x.runends,1)
  end
  return(item)
end

function pushfirst!(x::RLEVector{T1,T2},item) where {T1,T2 <: Integer}
  item = convert(T1,item) # Copying how base does it for arrays
  x.runends[:] = x.runends .+ 1
  if item != x.runvalues[1]
    pushfirst!(x.runvalues,item)
    pushfirst!(x.runends,1)
  end
  return(x)
end

function deleterun!(x::RLEVector,i::Integer)
    x.runends[i:end] = x.runends[i:end] .- widths(x,i)
    if (i > 1 && i < nrun(x) && x.runvalues[i-1] == x.runvalues[i+1])
        deleteat!(x.runvalues,(i-1):i)
        deleteat!(x.runends,(i-1):i)
    else
        deleteat!(x.runvalues,i)
        deleteat!(x.runends,i)
    end
    x
end

function decrement_run!(x::RLEVector,run::Integer)
  if widths(x,run) == 1
    deleterun!(x,run)
  else
    x.runends[run:end] .-= 1
  end
  return(x)
end

function deleteat!(x::RLEVector,i::Integer)
  run = ind2run(x,i)
  decrement_run!(x,run)
end

function insert!(x::RLEVector{T1,T2}, i::Integer, item) where {T1,T2 <: Integer}
    if i == 1
        pushfirst!(x,item)
    elseif i == length(x) + 1
        push!(x,item)
    else
        _item = convert(T1, item)
        (run, index_in_run, run_remainder) = ind2runcontext(x,i)
        if x.runvalues[run] == _item
            x.runends[run:end] = x.runends[run:end] .+ 1
        else
            new_nrun = nrun(x) + 2
            resize!(x.runvalues, new_nrun)
            resize!(x.runends, new_nrun)
            x.runvalues[(run + 2):end] = x.runvalues[run:(end-2)]
            x.runvalues[run + 1] = _item
            x.runends[(run + 2):end] = x.runends[run:(end-2)] .+ 1
            x.runends[run + 1] = i
            x.runends[run] = i - 1
        end
    end
    x
end

_default_splice = RLEVector(Union{}[],Int64[])
function splice!(x::RLEVector, i::Integer, ins::RLEVector=_default_splice)
    nrun_x = nrun(x)
    (1 <= i <= length(x)) || throw(BoundsError())
    if length(ins) == 0
        run = ind2run(x,i)
        current = x.runvalues[run]
        decrement_run!(x,run)
    else
        (run, index_in_run, run_remainder) = ind2runcontext(x,i)
        current = x.runvalues[run]
        # Splice ins into adjusted x
        widths!(x.runends)
        x.runends[run] = run_remainder
        nrun_out = run + nrun(ins) + ((nrun_x - run) + 1)
        resize!(x.runvalues, nrun_out)
        resize!(x.runends, nrun_out)
        x.runvalues[(run + (nrun_out - nrun_x)):nrun_out] = x.runvalues[run:nrun_x]
        x.runends[(run + (nrun_out - nrun_x)):nrun_out] = x.runends[run:nrun_x]
        x.runvalues[(run + 1):(run + nrun(ins))] = ins.runvalues
        x.runends[(run + 1):(run + nrun(ins))] = widths(ins)
        x.runends[run] = index_in_run - 1
        cumsum!(x.runends, x.runends)
        ree!(x.runvalues,x.runends)
    end
    return(current)
end

function splice!(x::RLEVector, index::UnitRange, ins::RLEVector=_default_splice)
    nrun_x = nrun(x)
    i_left = first(index)
    i_right = last(index)
    if i_left == i_right
        return(splice!(x,i_left,ins))
    end
    (run_left, run_right, index_in_run_left, run_remainder_right) = ind2runcontext(x,index)
    run_range = run_left:run_right
    if i_left > i_right # Insert without removing
        current = similar(x,0)
    else
        current = RLEVector(x.runvalues[run_range], x.runends[run_range] .- (i_left - 1))
        current.runends[end] = current.runends[end] - run_remainder_right
    end
    # Splice ins into adjusted x
    # FIXME: factor out code shared with setindex!(rle, rle, unit_range)
    widths!(x.runends)
    x.runends[run_right] = run_remainder_right
    nrun_out = run_left + nrun(ins) + ((nrun_x - run_right) + 1)
    resize!(x.runvalues, nrun_out)
    resize!(x.runends, nrun_out)
    x.runvalues[(run_right + (nrun_out - nrun_x)):nrun_out] = x.runvalues[run_right:nrun_x]
    x.runends[(run_right + (nrun_out - nrun_x)):nrun_out] = x.runends[run_right:nrun_x]
    x.runvalues[(run_left + 1):(run_left + nrun(ins))] = ins.runvalues
    x.runends[(run_left + 1):(run_left + nrun(ins))] = widths(ins)
    x.runends[run_left] = index_in_run_left - 1
    cumsum!(x.runends, x.runends)
    ree!(x.runvalues,x.runends)
    return(current)
end

function splice!(x::RLEVector, i::Union{Integer,UnitRange}, ins::AbstractArray)
  splice!(x,i,RLEVector(ins))
end

# Appended space initialized with zero unlike base array
function resize!(x::RLEVector, nl::Integer) # Based on base version for array
    l = length(x)
    if nl > l
        push!(x.runends,nl)
        push!(x.runvalues,0)
    else
        nl < 0 && throw(ArgumentError("new length must be ≥ 0"))
        (run, index_in_run, run_remainder) = ind2runcontext(x,nl)
        resize!(x.runends,run)
        resize!(x.runvalues,run)
        x.runends[end] = x.runends[end] - run_remainder
    end
    return(x)
end

function empty!(x::RLEVector)
  empty!(x.runvalues)
  empty!(x.runends)
  return(x)
end

### Set API
# Mostly handled by convert and promote rules in RLEVector-type.jl

function intersect(x::RLEVector, sets...)
    ok = trues(length(x.runvalues))
    for (i, v) in enumerate(x.runvalues)
        for s in sets
            if !in(v, s)
                ok[i] = false
                break
            end
        end
    end
    RLEVector( x.runvalues[ ok ], cumsum( widths(x)[ ok ] ) )
end


text:
function nilpotents(n)
    divs, _ = factor_m(n)
    multip = 1
    for x in divs
        multip *= x
    end
    return multip:multip:(n-1)
end


text:
using Logging

name_parts(loggername) = split(loggername, '.')

function configuration(configfile::String; section::String="")
    config = Pkg.TOML.parsefile(configfile)
    isempty(section) ? config : get_section(config, [section])
end

function get_section(config::Dict{String,Any}, path::Vector)
    isempty(path) && (return config)

    top = get(config, first(path)) do
        Dict{String,Any}()
    end

    (length(path) > 1) ? get_section(top, path[2:end]) : top
end

function flatten(config::Dict{String,Any}, path::Vector)
    result = Dict{String,Any}()
    for depth in 1:length(path)
        merge!(result, get_section(config, path[1:depth]))
    end
    result
end

function get_type(s::String, topmodule::Module=@__MODULE__)
    try
        T = topmodule
        for t in split(s, ".")
            T = Base.eval(T, Symbol(t))
        end
        return T
    catch ex
        nextmodule = parentmodule(topmodule)
        if nextmodule == topmodule
            error("Could not resolve logger $s. Ensure that the required logging packages are imported!")
        else
            return get_type(s, nextmodule)
        end
    end
end

function get_type(s::String, modules::Vector{Module})
    L = length(modules)
    for idx in 1:L
        try
            return get_type(s, modules[idx])
        catch ex
            (idx == L) && rethrow(ex)
        end
    end
end

logger(configfile::String, loggername::String; section::String="") = logger(configuration(configfile; section=section), name_parts(loggername))
logger(config::Dict{String,Any}, loggername::String) = logger(config, name_parts(loggername))
function logger(config::Dict{String,Any}, loggername::Vector)
    loggercfg = flatten(config, loggername)
    loggertypestr = loggercfg["type"]
    loggertopmodule = Base.eval(Main, Symbol(get(loggercfg, "topmodule", "Main")))
    modules = collect(Set([Main, loggertopmodule, @__MODULE__]))
    loggertype = get_type(loggertypestr, modules)

    logcompose(loggertype, config, loggercfg)
end

logcompose(::Type{T}, config::Dict{String,Any}, logger_config::Dict{String,Any}) where {T} = error("logcompose not implemented for type $T")

log_min_level(logger_config::Dict{String,Any}, default::String="Info") = getproperty(Logging, Symbol(get(logger_config, "min_level", default)))


text:
using LinearAlgebra, Plots

#-------------------------------------------------------------------------------------------
# FUNCTIONS

# Fresnel Coefficients for both TE and TM polarized light
function De(n, theta)
    return [ 1 1 ; n*cos(theta) -n*cos(theta)]
end
function Dm(n, theta)
    return [cos(theta) -cos(theta) ; n n ]
end

# determine k for each wavelength of visible light, incident on the surface at angle theta
function kFunc(n, λ, theta)
    """
    n: between ~1-10
    lambda: in nanometers (400-700)
    theta: angle for particular medium (in degrees)
    """
    # c0 = 2.998e+17;                     # nm/s
    return (n*2*pi)/(λ) * cos(theta)
end

# Propagation matrix
function P(n,d,λ,theta)
    k = kFunc(n,λ,theta)
    return [ exp(k*d*im) 0 ; 0 exp(-k*d*im)]
end

# stack the mirror
function mirror_stack(n1,n2,d1,d2,dSub,totalLayers,Polarization,λ)

    # Snell's law to determine angles in a repeating media
    theta1 = asin(sin(45)/n1);
    theta2 = asin(n1*sin(theta1)/n2);

    # Check polarization
    if Polarization == "TE"
        D = De
    elseif Polarization == "TM"
        D = Dm
    else
        println("Polarization must be either \"TE\" or \"TM\"")
    end

    # be very careful with number of layers here, can break
    # case 1 should have an even number of layers
    # case 2 should have an odd number of layers
    if totalLayers//2 == 0
        nSub = n1
        thetaSub = theta1
        return inv(D(nSub, dSub))*D(1, 45) * P(nSub, dSub, λ, thetaSub) * ( inv(D(n1,d1)) * D(n2,d2) * P(n2, d2, λ, theta2) * inv(D(n2,d2)) * D(n1,d1) * P(n1, d1, λ, theta1) )^((totalLayers)/2) * inv(D(1,45))*D(n1,theta1)
    else
        nSub = n2
        thetaSub = theta2 # inv(D(nSub, dSub))*D(1, 45) *
        return inv(D(nSub, dSub))*D(1, 45) * P(nSub, dSub, λ, thetaSub) * inv(D(n2,d2)) * D(n1,d1) * P(n1, d1, λ, theta1) * ( inv(D(n1,d1)) * D(n2,d2) * P(n2, d2, λ, theta2) * inv(D(n2,d2)) * D(n1,d1) * P(n1, d1, λ, theta1) )^((totalLayers - 1)/2) * inv(D(1,45))*D(n1,theta1)
    end
    # Consider: break the return lines into like 3 functions so it isn't so long
end

# determine the reflectivity of the light hitting the mirror
function reflect(n1,n2,d1,d2,dSub,totalLayers,Polarization,λ)
    matr = mirror_stack(n1,n2,d1,d2,dSub,totalLayers,Polarization,λ)
    return ((norm(matr[3]/matr[4])))^2
end

# plot it
function r_plot(n1,n2,d1,d2,dSub,totalLayers,Polarization,λ)
    # reflectivity across spectrum
    y = reflect.(n1, n2, d1, d2, dSub, totalLayers, Polarization, λ)

    # plot it
    plot(λ, y, ylims=(-0.1, 1.1), title="Reflectivity of $Polarization Polarized Light", xlabel="Wavelength [nm]", ylabel="Intensity [%]")

    savefig("reflectivity_" * "$Polarization" * ".png")
end
#-------------------------------------------------------------------------------------------
# TEST

# nm wavelength over visible spectrum
λ = 400:0.1:700;

# reflectivity
r_plot(1.46, 4.6, 52, 52, 1000, 21, "TM", λ);


text:
using StaticArrays
using Base: @_inline_meta
export reconstruct, DelayEmbedding, AbstractEmbedding, MTDelayEmbedding, embed
export WeightedDelayEmbedding

#####################################################################################
#                        Delay Embedding Reconstruction                             #
#####################################################################################
"""
    AbstractEmbedding
Super-type of embedding methods. Use `subtypes(AbstractEmbedding)` for available
methods.
"""
abstract type AbstractEmbedding <: Function end

"""
    DelayEmbedding(γ, τ) -> `embedding`
Return a delay coordinates embedding structure to be used as a functor,
given a timeseries and some index. Calling
```julia
embedding(s, n)
```
will create the `n`-th reconstructed vector of the embedded space, which has `γ`
temporal neighbors with delay(s) `τ`. See [`reconstruct`](@ref) for more.

**Be very careful when choosing `n`, because `@inbounds` is used internally.**
"""
struct DelayEmbedding{γ} <: AbstractEmbedding
    delays::SVector{γ, Int}
end

@inline DelayEmbedding(γ, τ) = DelayEmbedding(Val{γ}(), τ)
@inline function DelayEmbedding(::Val{γ}, τ::Int) where {γ}
    idxs = [k*τ for k in 1:γ]
    return DelayEmbedding{γ}(SVector{γ, Int}(idxs...))
end
@inline function DelayEmbedding(::Val{γ}, τ::AbstractVector) where {γ}
    γ != length(τ) && throw(ArgumentError(
    "Delay time vector length must equal the number of temporal neighbors."
    ))
    return DelayEmbedding{γ}(SVector{γ, Int}(τ...))
end

@generated function (r::DelayEmbedding{γ})(s::AbstractArray{T}, i) where {γ, T}
    gens = [:(s[i + r.delays[$k]]) for k=1:γ]
    quote
        @_inline_meta
        @inbounds return SVector{$γ+1,T}(s[i], $(gens...))
    end
end

# Weighted version
export WeightedDelayEmbedding
"""
    WeightedDelayEmbedding(γ, τ, w) -> `embedding`
Similar with [`DelayEmbedding`](@ref), but the entries of the
embedded vector are further weighted with `w^γ`.
See [`reconstruct`](@ref) for more.

**Be very careful when choosing `n`, because `@inbounds` is used internally.**
"""
struct WeightedDelayEmbedding{γ, T<:Real} <: AbstractEmbedding
    delays::SVector{γ, Int}
    w::T
end

@inline WeightedDelayEmbedding(γ, τ, w) = WeightedDelayEmbedding(Val{γ}(), τ, w)
@inline function WeightedDelayEmbedding(::Val{γ}, τ::Int, w::T) where {γ, T}
    idxs = [k*τ for k in 1:γ]
    return WeightedDelayEmbedding{γ, T}(SVector{γ, Int}(idxs...), w)
end

@generated function (r::WeightedDelayEmbedding{γ, T})(s::AbstractArray{X}, i) where {γ, T, X}
    gens = [:(r.w^($k) * s[i + r.delays[$k]]) for k=1:γ]
    quote
        @_inline_meta
        @inbounds return SVector{$γ+1,X}(s[i], $(gens...))
    end
end

"""
    reconstruct(s, γ, τ [, w])
Reconstruct `s` using the delay coordinates embedding with `γ` temporal neighbors
and delay `τ` and return the result as a [`Dataset`](@ref).

Use [`embed`](@ref) for the version that accepts the embedding dimension `D = γ+1`
instead.

## Description
### Single Timeseries
If `τ` is an integer, then the ``n``-th entry of the embedded space is
```math
(s(n), s(n+\\tau), s(n+2\\tau), \\dots, s(n+γ\\tau))
```
If instead `τ` is a vector of integers, so that `length(τ) == γ`,
then the ``n``-th entry is
```math
(s(n), s(n+\\tau[1]), s(n+\\tau[2]), \\dots, s(n+\\tau[γ]))
```

The reconstructed dataset can have same
invariant quantities (like e.g. lyapunov exponents) with the original system
that the timeseries were recorded from, for proper `γ` and `τ`.
This is known as the Takens embedding theorem [1, 2].
The case of different delay times allows reconstructing systems with many time scales,
see [3].

*Notice* - The dimension of the returned dataset (i.e. embedding dimension) is `γ+1`!

If `w` (a "weight") is provided as an extra argument, then the entries
of the embedded vector are further weighted with ``w^\\gamma``, like so
```math
(s(n), w*s(n+\\tau), w^2*s(n+2\\tau), \\dots,w^\\gamma * s(n+γ\\tau))
```

### Multiple Timeseries
To make a reconstruction out of a multiple timeseries (i.e. trajectory) the number
of timeseries must be known by type, so `s` can be either:

* `s::AbstractDataset{B}`
* `s::SizedAray{A, B}`

If the trajectory is for example ``(x, y)`` and `τ` is integer, then the ``n``-th
entry of the embedded space is
```math
(x(n), y(n), x(n+\\tau), y(n+\\tau), \\dots, x(n+γ\\tau), y(n+γ\\tau))
```
If `τ` is an `AbstractMatrix{Int}`, so that `size(τ) == (γ, B)`,
then we have
```math
(x(n), y(n), x(n+\\tau[1, 1]), y(n+\\tau[1, 2]), \\dots, x(n+\\tau[γ, 1]), y(n+\\tau[γ, 2]))
```

*Notice* - The dimension of the returned dataset is `(γ+1)*B`!

## References
[1] : F. Takens, *Detecting Strange Attractors in Turbulence — Dynamical
Systems and Turbulence*, Lecture Notes in Mathematics **366**, Springer (1981)

[2] : T. Sauer *et al.*, J. Stat. Phys. **65**, pp 579 (1991)

[3] : K. Judd & A. Mees, [Physica D **120**, pp 273 (1998)](https://www.sciencedirect.com/science/article/pii/S0167278997001188)
"""
function reconstruct(s::AbstractVector{T}, γ, τ) where {T}
    if γ == 0
        return Dataset{1, T}(s)
    end
    de::DelayEmbedding{γ} = DelayEmbedding(Val{γ}(), τ)
    return reconstruct(s, de)
end
function reconstruct(s::AbstractVector{T}, γ, τ, w) where {T}
    if γ == 0
        return Dataset{1, T}(s)
    end
    de = WeightedDelayEmbedding(Val{γ}(), τ, w)
    return reconstruct(s, de)
end
@inline function reconstruct(s::AbstractVector{T},
    de::Union{WeightedDelayEmbedding{γ}, DelayEmbedding{γ}}) where {T, γ}
    L = length(s) - maximum(de.delays)
    data = Vector{SVector{γ+1, T}}(undef, L)
    @inbounds for i in 1:L
        data[i] = de(s, i)
    end
    return Dataset{γ+1, T}(data)
end

"""
    embed(s, D, τ)
Perform a delay coordinates embedding on signal `s` with embedding dimension `D`
and delay time `τ`. The result is returned as a [`Dataset`](@ref), which is a
vector of static vectors.

See [`reconstruct`](@ref) for an advanced version that supports multiple delay
times and can reconstruct multiple timeseries efficiently.
"""
embed(s, D, τ) = reconstruct(s, D-1, τ)

#####################################################################################
#                              MultiDimensional R                                   #
#####################################################################################
"""
    MTDelayEmbedding(γ, τ, B) -> `embedding`
Return a delay coordinates embedding structure to be used as a functor,
given multiple timeseries (`B` in total), either as a [`Dataset`](@ref) or a
`SizedArray`), and some index.
Calling
```julia
embedding(s, n)
```
will create the `n`-th reconstructed vector of the embedded space, which has `γ`
temporal neighbors with delay(s) `τ`. See [`reconstruct`](@ref) for more.

**Be very careful when choosing `n`, because `@inbounds` is used internally.**
"""
struct MTDelayEmbedding{γ, B, X} <: AbstractEmbedding
    delays::SMatrix{γ, B, Int, X} # X = γ*B = total dimension number
end

@inline MTDelayEmbedding(γ, τ, B) = MTDelayEmbedding(Val{γ}(), τ, Val{B}())
@inline function MTDelayEmbedding(::Val{γ}, τ::Int, ::Val{B}) where {γ, B}
    X = γ*B
    idxs = SMatrix{γ,B,Int,X}([k*τ for k in 1:γ, j in 1:B])
    return MTDelayEmbedding{γ, B, X}(idxs)
end
@inline function MTDelayEmbedding(
    ::Val{γ}, τ::AbstractMatrix{<:Integer}, ::Val{B}) where {γ, B}
    X = γ*B
    γ != size(τ)[1] && throw(ArgumentError(
    "`size(τ)[1]` must equal the number of spatial neighbors."
    ))
    B != size(τ)[2] && throw(ArgumentError(
    "`size(τ)[2]` must equal the number of timeseries."
    ))
    return MTDelayEmbedding{γ, B, X}(SMatrix{γ, B, Int, X}(τ))
end
function MTDelayEmbedding(
    ::Val{γ}, τ::AbstractVector{<:Integer}, ::Val{B}) where {γ, B}
    error("Does not work with vector τ, only matrix or integer!")
end

@generated function (r::MTDelayEmbedding{γ, B, X})(
    s::Union{AbstractDataset{B, T}, SizedArray{Tuple{A, B}, T, 2, M}},
    i) where {γ, A, B, T, M, X}
    gensprev = [:(s[i, $d]) for d=1:B]
    gens = [:(s[i + r.delays[$k, $d], $d]) for k=1:γ for d=1:B]
    quote
        @_inline_meta
        @inbounds return SVector{$(γ+1)*$B,T}($(gensprev...), $(gens...))
    end
end

@inline function reconstruct(
    s::Union{AbstractDataset{B, T}, SizedArray{Tuple{A, B}, T, 2, M}},
    γ, τ) where {A, B, T, M}

    if γ == 0
        return Dataset{B, T}(s)
    end
    de::MTDelayEmbedding{γ, B, γ*B} = MTDelayEmbedding(γ, τ, B)
    reconstruct(s, de)
end
@inline function reconstruct(
    s::Union{AbstractDataset{B, T}, SizedArray{Tuple{A, B}, T, 2, M}},
    de::MTDelayEmbedding{γ, B, F}) where {A, B, T, M, γ, F}

    L = size(s)[1] - maximum(de.delays)
    X = (γ+1)*B
    data = Vector{SVector{X, T}}(undef, L)
    @inbounds for i in 1:L
        data[i] = de(s, i)
    end
    return Dataset{X, T}(data)
end

reconstruct(s::AbstractMatrix, args...) = reconstruct(Dataset(s), args...)


text:
execute = isempty(ARGS) || ARGS[1] == "run"

using Documenter
using Literate

# https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@example-block
ENV["GKSwstype"] = "100"
ENV["GKS_ENCODING"] = "utf-8"

# generate examples using Literate
lit = joinpath(@__DIR__, "lit")
src = joinpath(@__DIR__, "src")
gen = joinpath(@__DIR__, "src/generated")

base = "JuliaImageRecon/Examples"
repo_root_url =
    "https://github.com/$base/blob/main/docs/lit"
nbviewer_root_url =
    "https://nbviewer.org/github/$base/tree/gh-pages/generated"
binder_root_url =
    "https://mybinder.org/v2/gh/$base/gh-pages?filepath=generated"

for (root, _, files) in walkdir(lit), file in files
    splitext(file)[2] == ".jl" || continue # process .jl files only
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit => gen))[1]
    Literate.markdown(ipath, opath, documenter = execute; # run examples
        repo_root_url, nbviewer_root_url, binder_root_url)
    Literate.notebook(ipath, opath; execute = false, # no-run notebooks
        repo_root_url, nbviewer_root_url, binder_root_url)
end

# Documentation structure
ismd(f) = splitext(f)[2] == ".md"
pages(folder) =
    [joinpath("generated/", folder, f) for f in readdir(joinpath(gen, folder)) if ismd(f)]

isci = get(ENV, "CI", nothing) == "true"

format = Documenter.HTML(;
    prettyurls = isci,
    edit_link = "main",
    canonical = "https://juliaimagerecon.github.io/Examples",
#   assets = String[],
)

makedocs(;
    modules = Module[],
    authors = "Jeff Fessler and contributors",
    sitename = "Examples",
    format,
    pages = [
        "Home" => "index.md",
        "MRI" => pages("mri"),
        "CT" => pages("ct"),
    ],
)

if isci
    deploydocs(;
        repo = "github.com/JuliaImageRecon/Examples.git",
        devbranch = "main",
        versions = nothing,
        forcepush = true,
        push_preview = true,
    )
else
    @warn "may need to: rm -r src/generated/"
end

max_stars_repo_name: JuliaImageRecon/Examples

max_issues_repo_name: JuliaImageRecon/Examples

max_forks_repo_name: JuliaImageRecon/Examples


text:
function DX(delta, Nx, Ny)
    # This function creates a finite difference operator which represents a
    # first order x-derivitive.
    # INPUT
    # delta: resolution of the operator - scalar / a.u.
    # Nx: Number of elements in x-direction - integer / #
    # Ny: Number of elements in y-direction - integer / #
    # OUTPUT
    # operator:  (Nx*Ny X Nx*Ny)-Dim. FD-Dx operator - array / a.u.
    N = Nx*Ny
    next_vec = fill(1.0/(2.0*delta), N-1)
    if Ny > 1
        for ii in Nx:Nx:N-1 next_vec[ii] = 0.0 end
    end
    next_sparse_mat = sparse(2:N, 1:N-1, next_vec, N, N)
    operator = next_sparse_mat' - next_sparse_mat
    return operator
end

function DXX(delta, Nx, Ny)
    # This function creates a finite difference operator which represents a
    # second order x-derivitive.
    # INPUT
    # delta: resolution of the operator - scalar / a.u.
    # Nx: Number of elements in x-direction - integer / #
    # Ny: Number of elements in y-direction - integer / #
    # OUTPUT
    # operator:  (Nx*Ny X Nx*Ny)-Dim. FD-Dxx operator - array / a.u.
    N = Nx*Ny
    main_vec = fill(-2.0 /delta^2, N)
    next_vec = fill(1.0/delta^2, N-1)
    if Ny > 1
        for ii in Nx:Nx:(N-1) next_vec[ii] = 0.0 end
    end
    main_sparse_mat = sparse(1:N, 1:N, main_vec, N, N)
    next_sparse_mat = sparse(2:N, 1:(N-1), next_vec, N, N)
    operator = next_sparse_mat' + next_sparse_mat + main_sparse_mat
    return operator
end

function DY(delta, Nx, Ny)
    # This function creates a finite difference operator which represents a
    # first order y-derivitive.
    # INPUT
    # delta: resolution of the operator - scalar / a.u.
    # Nx: Number of elements in x-direction - integer / #
    # Ny: Number of elements in y-direction - integer / #
    # OUTPUT
    # operator:  (Nx*Ny X Nx*Ny)-Dim. FD-Dy operator - array / a.u.
    N = Nx*Ny
    next_vec = fill(1.0/(2*delta), (N-Ny))
    next_sparse_mat = sparse((1+Ny):N, 1:(N-Ny), next_vec, N, N)
    operator = next_sparse_mat' - next_sparse_mat
    return operator
end

function DYY(delta, Nx, Ny)
    # This function creates a finite difference operator which represents a
    # second order y-derivitive.
    # INPUT
    # delta: resolution of the operator - scalar / a.u.
    # Nx: Number of elements in x-direction - integer / #
    # Ny: Number of elements in y-direction - integer / #
    # OUTPUT
    # operator:  (Nx*Ny X Nx*Ny)-Dim. FD-Dyy operator - array / a.u.
    N = Nx*Ny
    main_vec = fill(-2.0/delta^2, N)
    next_vec = fill(1.0/delta^2, (N-Ny))
    main_sparse_mat = sparse(1:N, 1:N, main_vec, N, N)
    next_sparse_mat = sparse((1+Ny):N, 1:(N-Ny), next_vec, N, N)
    operator = next_sparse_mat + next_sparse_mat' + main_sparse_mat
    return operator
end


text:
function areaBC(r)
    return pi.*r.^2    # o ponto serve para calcular arrays
end


text:
#!/usr/bin/env julia

function kron_δ(i,j)
    if i == j
        return 1
    else
        return 0
    end
end

function harmean(i,j)
    2 / ((1/i)+(1/j))
end

function inv_degree(g,i)
    1/degree(g,i)
end

function strength(g,i)
    neigh = neighbors(g,i)
    strengthi = sum(map(x-> g.weights[x,i],neigh))
    strengtho = sum(map(x-> g.weights[i,x],neigh))
    strength = strengthi + strengtho
    strength
end

function strengthin(g,i)
    neigh = inneighbors(g,i)
    strength = sum(map(x-> g.weights[x,i],neigh))
    strength
end

function strengthout(g,i)
    neigh = outneighbors(g,i)
    strength = sum(map(x-> g.weights[x,i],neigh))
    strength
end

function inv_degree_ihara(g,i,v)
    try
        1/strength(g,i,v)
    catch
        0
    end
end

function rel_single(g,i)
    if outdegree(g,i) == 1
        1
    else
        0
    end
end

function mglobal_clustering_coefficient(g::AbstractGraph)
    c = 0
    ntriangles = 0
    for v in vertices(g)
        neighs = neighbors(g, v)
        for i in neighs, j in neighs
            i == j && continue
            if has_edge(g, i, j)
                c += 1
            end
        end
        k = degree(g, v)
        ntriangles += k * (k - 1)
    end
    ntriangles == 0 && return 0.
    return c / ntriangles
end

function inv_degree_N(g,i)
    if degree(g,i) != 1
        1/(degree(g,i)-1)
    else
        1/degree(g,i)
    end
end

function nonBM_embedding(g)
  NB = la_NBM(g)
  λ, eigve = eigs(NB, nev=20)
  los_reales = Array{Float64}(0)
  los_indices_reales = Array{Int64}(0)
  for i in 1:length(λ)
        if imag(λ[i]) == 0 && real(λ[i]) != 0 && i*real(λ[i]) >= sqrt(real(λ[1]))
          push!(los_reales,λ[i])
          push!(los_indices_reales,i)
      end
  end
  matrix_embedded = sub(real(eigve),:,los_indices_reales)
  length(los_reales), los_indices_reales
  ϕ = zeros(Float64, nv(g), length(los_reales))
  for n=1:(length(los_reales))
      v= matrix_embedded[:,n]
      ϕ[:,n] = contract(g, v)
  end
  return length(los_reales), ϕ
end

function grados(g)
	grados = Array{Int64}(0)
	gradosdiv = Array{Int64}(0)
	for i in vertices(g)
		push!(grados,degree(g,i))
	end
	grados
end

function avg_grados(g)
    grado = grados(g)
    avg_grad = mean(grado)
    avg_grad
end

function max_grados(g)
    grado = grados(g)
    maxi = maximum(grado)
    maxi
end

function mat_an(g,v,NBM)
	unos = ones(2ne(g))/sqrt(2ne(g))
	unnos = unos*unos'
	R = NBM - unnos
	R
end

function eigens(M,threshold)
	valores, vectores = eigs(M,which=:LR,nev=3)
    ind = findall(x->imag(x)==0.0,valores)
    valores = valores[ind]
    vectores = vectores[:,ind]
    if real(last(valores)) > threshold
		valores, vectores = eigs(M,which=:LR,nev=4)
        ind = findall(x->imag(x)==0.0,valores)
        valores = valores[ind]
        vectores = vectores[:,ind]
		if real(last(valores)) > threshold
		    valores, vectores = eigs(M,which=:LR,nev=5)
            ind = findall(x->imag(x)==0.0,valores)
            valores = valores[ind]
            vectores = vectores[:,ind]
		    if real(last(valores)) > threshold
		        valores, vectores = eigs(M,which=:LR,nev=10)
                ind = findall(x->imag(x)==0.0,valores)
                valores = valores[ind]
                vectores = vectores[:,ind]
		        if real(last(valores)) > threshold
		            valores, vectores = eigs(M,which=:LR,nev=nv(g))
                    ind = findall(x->imag(x)==0.0,valores)
                    valores = valores[ind]
                    vectores = vectores[:,ind]
		        end
		    end
		end
	end
	valores, vectores
end

function num_com(v,threshold)
	cuantos = Array{Float64}(undef,0)
	index = Array{Int64}(undef,0)
	for i in 1:length(v)
		if imag(v[i]) == 0 && real(v[i]) > threshold
		    push!(cuantos,real(v[i]))
		    push!(index,i)
		end
	end
    cuantos,index
end

function embe_nbm(g,ind,vec)
	reales = real(vectores[:,index])
	up = reales[1:nv(g),:]
	down = reales[nv(g)+1:2*nv(g),:]
	matrix_embedded = up + down
	matrix_embedded
end

function contract(g,v,ei)
    y = zeros(Float64, nv(g))
    for i in 1:nv(g)
        for j in inneighbors(g,i)
            u = ei[Edge(j,i)]
            y[i] += v[u]
        end
    end
    y
end

function contraction(g,vec,emb,ei)
	contracted = zeros(Float64, nv(g), length(vec))
	for n in 1:length(vec)
    	contracted[:,n] = contract(g,emb,ei)
	end
	contracted
end