text:
function tcm_estimateForest_RTS_kcv_weighted(N, K_sample, d_list, nfolds, transaction_folds, transactions, time_limit_EM, perTol_EM, initial_left, initial_right, initial_product, initial_isLeaf)
	unique_fold_ids = unique(transaction_folds)

	KL_total_holdout_fold = zeros(length(d_list), nfolds)
	noinf_count_holdout_fold = zeros(length(d_list), nfolds)

	start_time = time()

	for fold_ind in 1:nfolds
		fold = unique_fold_ids[fold_ind]

		@show fold

		@show unique_fold_ids

		println("Pre-subsetting ")

		#@show transaction_folds
		#temp = find(transaction_folds .!= fold)
		#@show temp
		@show size(transaction_folds)
		@show size(transactions)
		@show fold

		temp = find(transaction_folds .!= fold)
		@show temp[1:5]
		trainfolds_transactions = transactions[ find(transaction_folds .!= fold), :]
		println("Post-subsetting")
		testfold_transactions = transactions[ find(transaction_folds .== fold), :]

		trainfolds_assortments, trainfolds_transaction_counts = tcm_transactionsToCounts(N, trainfolds_transactions)
		testfold_assortments, testfold_transaction_counts = tcm_transactionsToCounts(N, testfold_transactions)

		for i in 1:length(d_list)
			println("Depth ", d_list[i], ", fold ", fold)
			depth_limit = d_list[i]

			if (isempty(initial_left))
				lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time = tcm_estimateForest_RTS(N, trainfolds_assortments, trainfolds_transaction_counts, K_sample, depth_limit, perTol_EM, time_limit_EM)

			else
				lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time = tcm_estimateForest_RTS_withInitial(N, trainfolds_assortments, trainfolds_transaction_counts, K_sample, depth_limit, perTol_EM, time_limit_EM,
																																initial_left, initial_right, initial_product, initial_isLeaf)
			end

			K = length(lambda)
			testfold_forest_predict, testfold_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, testfold_assortments)

			testfold_KL, testfold_noinf_count, testfold_KL_avg = tcm_evaluateKL(N, testfold_forest_predict, testfold_transaction_counts)

			KL_total_holdout_fold[i, fold_ind] = testfold_KL #testfold_KL / sum(testfold_transaction_counts)
			noinf_count_holdout_fold[i, fold_ind] = testfold_noinf_count
		end
	end

	end_time = time()
	elapsed_time = end_time - start_time

	KL_avg_by_d = sum(KL_total_holdout_fold, 2) ./ sum(noinf_count_holdout_fold, 2)  * 100

	@show KL_total_holdout_fold
	@show noinf_count_holdout_fold
	@show KL_avg_by_d

	best_i = indmin(KL_avg_by_d)
	best_d = d_list[best_i]
	best_KL_avg_cv = KL_avg_by_d[best_i]

	return best_d, best_KL_avg_cv, elapsed_time

end

function tcm_estimateForest_RTS_kcv_weighted_ass(N, K_sample, d_list, nfolds, assortment_folds, assortments, transaction_counts, time_limit_EM, perTol_EM, initial_left, initial_right, initial_product, initial_isLeaf)
	unique_fold_ids = unique(assortment_folds)

	KL_total_holdout_fold = zeros(length(d_list), nfolds)
	noinf_count_holdout_fold = zeros(length(d_list), nfolds)

	start_time = time()

	for fold_ind in 1:nfolds
		fold = unique_fold_ids[fold_ind]

		# trainfolds_assortments, trainfolds_transaction_counts = tcm_transactionsToCounts(N, trainfolds_transactions)
		# testfold_assortments, testfold_transaction_counts = tcm_transactionsToCounts(N, testfold_transactions)

		trainfolds_assortments = assortments[ find(assortment_folds .!= fold), :]
		testfold_assortments = assortments[ find(assortment_folds .== fold), :]

		trainfolds_transaction_counts = transaction_counts[ find(assortment_folds .!= fold), :]
		testfold_transaction_counts = transaction_counts[ find(assortment_folds .== fold), :]

		for i in 1:length(d_list)
			println("Depth ", d_list[i], ", fold ", fold)
			depth_limit = d_list[i]

			if (isempty(initial_left))
				lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time = tcm_estimateForest_RTS(N, trainfolds_assortments, trainfolds_transaction_counts, K_sample, depth_limit, perTol_EM, time_limit_EM)

			else
				lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time = tcm_estimateForest_RTS_withInitial(N, trainfolds_assortments, trainfolds_transaction_counts, K_sample, depth_limit, perTol_EM, time_limit_EM,
																																initial_left, initial_right, initial_product, initial_isLeaf)
			end

			K = length(lambda)
			testfold_forest_predict, testfold_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, testfold_assortments)

			testfold_KL, testfold_noinf_count, testfold_KL_avg = tcm_evaluateKL(N, testfold_forest_predict, testfold_transaction_counts)

			KL_total_holdout_fold[i, fold_ind] = testfold_KL #testfold_KL / sum(testfold_transaction_counts)
			noinf_count_holdout_fold[i, fold_ind] = testfold_noinf_count
		end
	end

	end_time = time()
	elapsed_time = end_time - start_time

	KL_avg_by_d = sum(KL_total_holdout_fold, 2) ./ sum(noinf_count_holdout_fold, 2)  * 100

	best_i = indmin(KL_avg_by_d)
	best_d = d_list[best_i]
	best_KL_avg_cv = KL_avg_by_d[best_i]

	return best_d, best_KL_avg_cv, elapsed_time

end

function tcm_estimateForest_RTS_RMredo_kcv_weighted_ass(N, K_sample, d_list, nfolds, assortment_folds, assortments, transaction_counts, time_limit, perTol, perTol_EM)
	unique_fold_ids = unique(assortment_folds)

	KL_total_holdout_fold = zeros(length(d_list), nfolds)
	noinf_count_holdout_fold = zeros(length(d_list), nfolds)

	start_time = time()

	for fold_ind in 1:nfolds
		fold = unique_fold_ids[fold_ind]

		# trainfolds_assortments, trainfolds_transaction_counts = tcm_transactionsToCounts(N, trainfolds_transactions)
		# testfold_assortments, testfold_transaction_counts = tcm_transactionsToCounts(N, testfold_transactions)

		trainfolds_assortments = assortments[ find(assortment_folds .!= fold), :]
		testfold_assortments = assortments[ find(assortment_folds .== fold), :]

		trainfolds_transaction_counts = transaction_counts[ find(assortment_folds .!= fold), :]
		testfold_transaction_counts = transaction_counts[ find(assortment_folds .== fold), :]

		for i in 1:length(d_list)
			println("Depth ", d_list[i], ", fold ", fold)
			depth_limit = d_list[i]

			ranking_lambda, orderings, ranking_loglik, ranking_elapsed_time = tcm_estimateRanking(N, trainfolds_assortments, trainfolds_transaction_counts, time_limit, perTol)
			initial_left, initial_right, initial_product, initial_isLeaf = tcm_convertRankingToForest(N, orderings)

			lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time = tcm_estimateForest_RTS_withInitial(N, trainfolds_assortments, trainfolds_transaction_counts, K_sample, depth_limit, perTol_EM, time_limit,
																															initial_left, initial_right, initial_product, initial_isLeaf)

			K = length(lambda)
			testfold_forest_predict, testfold_forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, testfold_assortments)

			testfold_KL, testfold_noinf_count, testfold_KL_avg = tcm_evaluateKL(N, testfold_forest_predict, testfold_transaction_counts)

			KL_total_holdout_fold[i, fold_ind] = testfold_KL #testfold_KL / sum(testfold_transaction_counts)
			noinf_count_holdout_fold[i, fold_ind] = testfold_noinf_count
		end
	end

	end_time = time()
	elapsed_time = end_time - start_time

	KL_avg_by_d = sum(KL_total_holdout_fold, 2) ./ sum(noinf_count_holdout_fold, 2)  * 100

	best_i = indmin(KL_avg_by_d)
	best_d = d_list[best_i]
	best_KL_avg_cv = KL_avg_by_d[best_i]

	return best_d, best_KL_avg_cv, elapsed_time

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
@testset "args" begin
    s = """
        a "bc def ghi" kl 1
        """
    args = FP.split_args(s)
    args[1] // "a"
    args[2] // "\"bc def ghi\""
    args[3] // "kl"
    args[4] // "1"

    a = "abc def" |> FP.split_args
    a[1] // "abc"
    a[2] // "def"

    a = "" |> FP.split_args
    isempty(a)

    a = "abc" |> FP.split_args
    a[1] // "abc"

    a = "\"abc\"" |> FP.split_args
    a[1] // "\"abc\""

    a = "abc \"def\" \"\"\"ghi\"\"\" klm" |> FP.split_args
    a[1] // "abc"
    a[2] // "def"
    a[3] // "\"ghi\""
    a[4] // "klm"
end

ext: jl


text:
include("../src/EDKit.jl")
using .EDKit
using LinearAlgebra, Test

#-------------------------------------------------------------------------------------------------------------------------
# SO(3) Random Quasi-symmetric Model
#-------------------------------------------------------------------------------------------------------------------------
@testset "SO(3) qsymm translation" begin
    L = 6
    h = rand()
    rh = begin
        vecs = zeros(ComplexF64, 9, 4)
        vecs[8,1] = 1
        vecs[6,1] = -1
        vecs[7,2] = 1
        vecs[3,2] = -1
        vecs[4,3] = 1
        vecs[2,3] = -1
        vecs[3,4] = 1
        vecs[5,4] = -1
        vecs[7,4] = 1
        rm = rand(ComplexF64, 4,4) .- 0.5 |> Hermitian |> Array
        vecs * rm * vecs'
    end
    H  = trans_inv_operator(rh, 2, L)
    H += trans_inv_operator(spin((h, "z"), D=3), 1, L)
    E, V = Array(H) |> Hermitian |> eigen

    for inds in Any[1:L÷2, [1,2], [3,5], [2,4,6]]
        ts = Float64[]
        for i = 0:L-1
            tb = TranslationalBasis(i, L, base = 3)
            H  = trans_inv_operator(rh, 2, tb)
            H += trans_inv_operator(spin((h, "z"), D=3), 1, tb)
            e, v = Array(H) |> Hermitian |> eigen
            ee = [ent_S(v[:, i], inds, tb) for i=1:size(v, 2)]
            append!(ts, ee)
        end
        EE = [ent_S(V[:, i], inds, TensorBasis(L, base=3)) for i=1:size(V, 2)]
        @test norm(sort(EE) - sort(ts)) ≈ 0.0 atol = 1e-7
    end
end

@testset "SO(3) qsymm translation + inversion" begin
    L = 6
    h = rand()
    perm = [1, 4, 7, 2, 5, 8, 3, 6, 9]
    rh = begin
        vecs = zeros(ComplexF64, 9, 4)
        vecs[8,1] = 1
        vecs[6,1] = -1
        vecs[7,2] = 1
        vecs[3,2] = -1
        vecs[4,3] = 1
        vecs[2,3] = -1
        vecs[3,4] = 1
        vecs[5,4] = -1
        vecs[7,4] = 1
        rm = rand(ComplexF64, 4,4) .- 0.5 |> Hermitian |> Array
        mat = vecs * rm * vecs'
        mat + mat[perm, perm]
    end

    function solve(basis)
        H  = trans_inv_operator(rh, 2, basis)
        H += trans_inv_operator(spin((h, "z"), D=3), 1, basis)
        e, v = Array(H) |> Hermitian |> eigen
        v
    end
    EE(v, inds, basis) = [ent_S(v[:, i], inds, basis) for i=1:size(v, 2)]

    # K = 0
    ba = TranslationalBasis(0, L, base=3)
    be = TranslationParityBasis(0, +1, L, base=3)
    bo = TranslationParityBasis(0, -1, L, base=3)
    va, ve, vo = solve(ba), solve(be), solve(bo)
    for inds in Any[1:L÷2, [1,2], [3,5], [2,4,6]]
        eea = EE(va, inds, ba)
        eee = EE(ve, inds, be)
        eeo = EE(vo, inds, bo)
        @test norm(sort(eea) - sort(vcat(eee, eeo))) ≈ 0.0 atol = 1e-7
    end

    # K = π
    ba = TranslationalBasis(L÷2, L, base=3)
    be = TranslationParityBasis(L÷2, +1, L, base=3)
    bo = TranslationParityBasis(L÷2, -1, L, base=3)
    va, ve, vo = solve(ba), solve(be), solve(bo)
    for inds in Any[1:L÷2, [1,2], [3,5], [2,4,6]]
        eea = EE(va, inds, ba)
        eee = EE(ve, inds, be)
        eeo = EE(vo, inds, bo)
        @test norm(sort(eea) - sort(vcat(eee, eeo))) ≈ 0.0 atol = 1e-7
    end
end

#-------------------------------------------------------------------------------------------------------------------------
# Random model with spin-flip symmetry
#-------------------------------------------------------------------------------------------------------------------------
@testset "Random spin-flip" begin
    L = 6
    perm = [9,8,7,6,5,4,3,2,1]
    rh = begin
        mat = rand(ComplexF64, 9, 9) |> Hermitian |> Array
        mat + mat[perm, perm]
    end

    function solve(basis)
        H  = trans_inv_operator(rh, 2, basis)
        e, v = Array(H) |> Hermitian |> eigen
        v
    end
    EE(v, inds, basis) = [ent_S(v[:, i], inds, basis) for i=1:size(v, 2)]

    for i = 0:L-1
        ba = TranslationalBasis(i, L, base=3)
        be = TranslationFlipBasis(i, +1, L, base=3)
        bo = TranslationFlipBasis(i, -1, L, base=3)
        va, ve, vo = solve(ba), solve(be), solve(bo)
        for inds in Any[1:L÷2, [1,2], [3,5], [2,4,6]]
            eea = EE(va, inds, ba)
            eee = EE(ve, inds, be)
            eeo = EE(vo, inds, bo)
            @test norm(sort(eea) - sort(vcat(eee, eeo))) ≈ 0.0 atol = 1e-7
        end
    end

end


text:
function _worker_next_step(time::DateTime, person::Person)
    state = person.movement_state

    # Handle day end
    if Dates.hour(time) == DAY_END
        person.movement_state = 1
        return true
    end

    if state == 1 # Morning state
        person.movement_state = 2
    elseif state == 2 # Transport before work state
        person.movement_state = 3
    elseif state == 3 && _is_working_end_time(time) # Working state
        person.movement_state = 4
    elseif state == 4 # Transport after work state
        person.movement_state = 5
    elseif state == 5 # Home after work
        p = rand()
        if p < _worker_go_leisure(time)
            person.movement_state = 6
        end
    elseif state == 6 # Leisure place
        p = rand()
		person.time_spent_out += 1
        if p < _come_back_from_leisure(person.time_spent_out)
            person.movement_state = 5
			person.time_spent_out = 0
        end
    end

    return state != person.movement_state
end

function _student_next_step(time::DateTime, person::Person)
    state = person.movement_state

    # Handle day end
    if Dates.hour(time) == DAY_END
        person.movement_state = 1
        return true
    end

    if state == 1 # Morning state
        person.movement_state = 2
    elseif state == 2 # Transport before work state
        person.movement_state = 3
    elseif state == 3 && _is_school_end_time(time) # Working state
        person.movement_state = 4
    elseif state == 4 # Transport after work state
        person.movement_state = 5
    elseif state == 5 # Home after work
        p = rand()
        if p < _student_go_leisure(time)
            person.movement_state = 6
        end
    elseif state == 6 # Leisure place
        p = rand()
		person.time_spent_out += 1
        if p < _come_back_from_leisure(person.time_spent_out)
            person.movement_state = 5
			person.time_spent_out = 0
        end
    end

    return state != person.movement_state
end

function _home_next_step(time::DateTime, person::Person)
    state = person.movement_state

    # Handle day end
    if Dates.hour(time) == DAY_END
        person.movement_state = 1
        return true
    end

    p = rand()
    if state == 1
		if p < _home_go_leisure(time)
        	person.movement_state = 2
		end
    elseif state == 2
		person.time_spent_out += 1
        if p < _come_back_from_leisure(person.time_spent_out)
            person.movement_state = 1
            person.time_spent_out = 0
        end
    end

    return state != person.movement_state
end

function _worker_go_leisure(time::DateTime)
    return (WORK_END + 1)^5 * 0.3 / Dates.hour(time)^5
end

function _student_go_leisure(time::DateTime)
    return pdf(Normal(20, 1.5), Dates.hour(time)) + 0.3
end

"""
	Get probability go to leisure for home persons, based on the day phase (morning, afternoon, evening)
"""
function _home_go_leisure(time::DateTime)
	μ = 0.0
	σ = 0.0
	if _is_morning(time)
		μ = mean([DAY_START, DAY_MEDIAN])
		σ = 0.8
	elseif _is_afternoon(time)
		μ = mean([DAY_MEDIAN, WORK_END])
		σ = 0.8
	elseif _is_evening(time)
		μ = mean([WORK_END, DAY_END])
		σ = 0.8
	end

	return pdf(Normal(μ, σ), Dates.hour(time))
end

"""
	Get the probability of stay in the current leisure place after x hour
	Made to have the following state:
		f(1) ≈ 0.4
		f(2) ≈ 0.2
		f(3) ≈ 0.1
		f(4) ≈ 0
"""
function _stay_leisure(x::Int)
    return -0.3153951 + (1132000 - -0.3153951) / (1 + (x / 1.252143e-7)^0.8598911)
end

function _come_back_from_leisure(time_out::Int)
    return round(1 - _stay_leisure(time_out), digits = 2)
end

# function _worker_get_probability_stay_leisure(time::DateTime)
# 	return (WORK_END + 2)^10 * 0.8 / Dates.hour(time)^10
# end

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


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

using Markdown
using InteractiveUtils

# ╔═╡ 3415208a-ab79-11ec-1f49-bdfda5782527
begin
	using StatsPlots
	using LinearAlgebra
	using QuadGK
md"""
# HW 10

Eric Nguyen
"""
end

# ╔═╡ e9dba4a1-2b85-4f5d-90fc-87ff6a4317ab
md"""
### Problem 1
"""

# ╔═╡ bb4bf5a6-c482-42bc-8b4c-e8ae7b48bb1d
md"""
$\underline{\underline{A}} = \begin{bmatrix} 1 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 0 & -1 \end{bmatrix}$

We write the system of equations ``\underline{\underline{A}} \, \underline{x} = \underline{c}`` as

$\begin{align*}
x_1 + x_2 + x_3 &= c_1 \\
x_1 + x_3 &= c_2 \\
x_1 - x_3 &= c_3
\end{align*}$

From this system of equations, we see that

$x_1 = x_3 + c_3 \implies x_3 = \frac{c_2 - c_3}{2} \implies x_1 = c_3 + \frac{c_2 - c_3}{2} = \frac{c_2 + c_3}{2}$

and

$x_2 = c_1 - \frac{c_2 - c_3}{2} - \frac{c_2 + c_3}{2} = c_1 - c_2$

so

$\underline{x} = \begin{bmatrix} \frac{1}{2} (c_2 + c_3) \\ c_1 - c_2 \\ \frac{1}{2} (c_2 - c_3) \end{bmatrix} = \begin{bmatrix} \frac{1}{2} c_2 + \frac{1}{2} c_3 \\ c_1 - c_2 \\ \frac{1}{2} c_2 - \frac{1}{2} c_3 \end{bmatrix} = \begin{bmatrix} 0 & \frac{1}{2} & \frac{1}{2} \\ 1 & -1 & 0 \\ 0 & \frac{1}{2} & -\frac{1}{2} \end{bmatrix} \underline{c} = \underline{\underline{A}}^{-1} \underline{c}$

which implies that

$\underline{\underline{A}}^{-1} = \begin{bmatrix} 0 & \frac{1}{2} & \frac{1}{2} \\ 1 & -1 & 0 \\ 0 & \frac{1}{2} & -\frac{1}{2} \end{bmatrix}$
"""

# ╔═╡ c9a48dd7-2e94-4e44-80ff-a0b7966f5a2f
md"""
### Problem 2

We want to approximate the function ``f(x) = e^x`` over the domain ``-1 ≤ x ≤ 1`` in the form

$f(x) = \sum_{ℓ = 0}^∞ A_ℓ P_ℓ(x)$

The Legendre Polynomials ``P_ℓ(x)`` obey the orthogonality relationship

$∫_{-1}^1 P_n(x) P_m(x) \,dx = \frac{2}{2n + 1} \delta_{nm}$

Integrating the product of the function and a Legendre polynomial from ``-1`` to ``1`` gives us

$∫_{-1}^1 f(x) P_ℓ(x) \,dx = ∫_{-1}^1 \left(\sum_{ℓ = 0}^∞ A_ℓ P_ℓ(x)\right) P_ℓ(x) \,dx$

$= A_ℓ ∫_{-1}^1 P_ℓ(x) P_ℓ(x) \,dx = A_ℓ \left(\frac{2}{2ℓ + 1} \delta_{ℓℓ}\right)$

$\implies A_ℓ = \frac{2ℓ + 1}{2} ∫_{-1}^1 f(x) P_ℓ(x) \,dx$

I will use the first four terms.
The first four Legendre polynomials are

$\begin{align*}
P_0(x) &= 1 \\
P_1(x) &= x \\
P_2(x) &= \frac{1}{2} (3x^2 - 1) \\
P_3(x) &= \frac{1}{2} (5x^3 - 3x)
\end{align*}$

So the first four values of ``A_ℓ`` are

$\begin{align*}
A_0 &= \frac{1}{2} ∫_{-1}^1 e^x \,dx = 1.1752 \\
A_1 &= \frac{3}{2} ∫_{-1}^1 e^x ⋅ (x) \,dx = 1.10364 \\
A_2 &= \frac{5}{2} ∫_{-1}^1 e^x ⋅ \left[\frac{1}{2} (3x^2 - 1)\right] \,dx = 0.357814 \\
A_3 &= \frac{7}{2} ∫_{-1}^1 e^x ⋅ \left[\frac{1}{2} (5x^3 - 3x)\right] \,dx = 0.0704556
\end{align*}$

and we represent the approximation of ``f(x) = e^x`` as

$f(x) = 1.1752 P_0(x) + 1.10364 P_1(x) + 0.357814 P_2(x) + 0.0704556 P_3(x)$
"""

# ╔═╡ 92098eba-7e9d-4d07-8e38-8599778e4034
let
	f(x) = exp(x)
	P0(x) = 1
	P1(x) = x
	P2(x) = (1/2) * (3x^2 - 1)
	P3(x) = (1/2) * (5x^3 - 3x)
	A(ℓ,P) = quadgk(x -> ((2ℓ + 1) / 2) * f(x) * P(x), -1, 1)
	A = [x for (x,y) ∈ [A(0,P0), A(1,P1), A(2,P2), A(3,P3)]]

	F(x) = A[1] * P0(x) + A[2] * P1(x) + A[3] * P2(x) + A[4] * P3(x)
	plot(exp, xlims=[-1,1], label="f(x) = eˣ", legend=:topleft, title="Approximation of f(x) = eˣ compared to exact function")
	plot!(F, label="Approximation")
end

# ╔═╡ bef51e5b-b663-4139-8a17-018b473ca188
md"""
The plot shows that the approximation is nearly identical to the exact function for ``-1 ≤ x ≤ 1``.
"""

# ╔═╡ 201e9fed-4e7b-480e-a523-e879ab34ec2f
md"""
### Problem 3

$\text{tr}(\underline{\underline{A}} \,\underline{\underline{B}}) = [\underline{\underline{A}} \,\underline{\underline{B}}]_{ii} = \sum_{k = 1}^N A_{ik} B_{ki} = \sum_{k = 1}^N B_{ik} A_{ki} = [\underline{\underline{B}} \,\underline{\underline{A}}]_{ii} = \text{tr}(\underline{\underline{B}} \,\underline{\underline{A}})$
"""

# ╔═╡ 50258060-a9b6-468c-a57a-c74c9444fe4a
md"""
### Problem 4

**(1).**
Rotations in the 3D space around the ``z``-axis are given by the following matrix,

$\begin{bmatrix} \cos{\phi} & -\sin{\phi} & 0 \\ \sin{\phi} & \cos{\phi} & 0 \\ 0 & 0 & 1 \end{bmatrix}$

So,

$\underline{\underline{A}} = \begin{bmatrix} \cos{90^\circ} & -\sin{90^\circ} & 0 \\ \sin{90^\circ} & \cos{90^\circ} & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}$

**(b).**
Rotations in the 3D space around the ``x``-axis are given by the following matrix,

$\begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos{\phi} & -\sin{\phi} \\ 0 & \sin{\phi} & \cos{\phi} \end{bmatrix}$

So,

$\underline{\underline{B}} =
\begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos{90^\circ} & -\sin{90^\circ} \\ 0 & \sin{90^\circ} & \cos{90^\circ} \end{bmatrix} =
\begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{bmatrix}$

**(c).**

$\underline{\underline{A}} \,\underline{\underline{B}} = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}$

$\underline{\underline{B}} \,\underline{\underline{A}} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & -1 & 0 \\ 0 & 0 & -1 \\ 1 & 0 & 0 \end{bmatrix}$

So,

$\underline{\underline{A}} \,\underline{\underline{B}} ≠ \underline{\underline{B}} \,\underline{\underline{A}}$

**(d).**

$\underline{\underline{A}} \,\underline{\underline{B}} \,\underline{v} = \begin{bmatrix} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}$

$\underline{\underline{B}} \,\underline{\underline{A}} \,\underline{v} = \begin{bmatrix} 0 & -1 & 0 \\ 0 & 0 & -1 \\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} = \begin{bmatrix} -1 \\ 0 \\ 0 \end{bmatrix}$
"""

# ╔═╡ 7492da11-24f6-4e04-84a5-831e615c4504
let
	ϕ = 90
	A = [cosd(ϕ) -sind(ϕ) 0; sind(ϕ) cosd(ϕ) 0; 0 0 1]
	B = [1 0 0; 0 cosd(ϕ) -sind(ϕ); 0 sind(ϕ) cosd(ϕ)]
	v = [0;1;0]
	O = (0,0,0)
	B*A*v
	scatter([O], label=false, xlims=[-1,1], ylims=[-1,1], zlims=[-1,1], xlabel="x-axis", ylabel="y-axis", zlabel="z-axis", color=:black)
	plot!([O, Tuple(A*B*v)], label="ABv", lw=2, color=:blue)
	plot!([O, Tuple(B*A*v)], label="BAv", lw=2, color=:red)
end

# ╔═╡ 36330a8c-e1c6-46ee-8255-b583e5cba6ff
md"""
### Problem 5

$\underline{v} = \begin{bmatrix} i \\ 1 \\ -1 \end{bmatrix} \qquad \underline{u} = \begin{bmatrix} 2i \\ -2 \\ 0 \end{bmatrix}$

The norm of ``\underline{v}`` is

$\sqrt{⟨v \vert v⟩} = \sqrt{v_i^* v_i} = \sqrt{(-i)(i) + (1)(1) + (-1)(-1)} = \sqrt{1 + 1 + 1} = \sqrt{3}$

The norm of ``\underline{u}`` is

$\sqrt{⟨u \vert u⟩} = \sqrt{u_i^* u_i} = \sqrt{(-2i)(2i) + (-2)(-2) + (0)(0)} = \sqrt{4 + 4 + 0} = \sqrt{8} = 2 \sqrt{2}$

We can show that ``\underline{v}`` and ``\underline{u}`` are orthogonal to each other by showing that their inner product equals zero.

$⟨v \vert u⟩ = \underline{\tilde{v}} \,\underline{u} = v_i^* u_i = (-i)(2i) + (1)(-2) + (-1)(0) = 2 - 2 = 0$

To find a vector ``\underline{w}`` with unit norm orthogonal to both ``\underline{v}`` and ``\underline{u}``, we can just take the cross product ``\underline{v} × \underline{u}`` and normalize it.

$\underline{w} = \underline{v} × \underline{u} = \begin{bmatrix} (1)(0) - (-2)(-1) \\ (2i)(-1) - (i)(0) \\ (i)(-2) - (2i)(1) \end{bmatrix} = \begin{bmatrix} -2 \\ -2i \\ -4i \end{bmatrix}$

We normalize ``\underline{w}`` by dividing it by its norm.
The norm of ``\underline{w}`` is

$\sqrt{⟨w,w⟩} = \sqrt{w_i^* w_i} = \sqrt{(-2)(-2) + (2i)(-2i) + (4i)(-4i)} = \sqrt{4 + 4 + 16} = \sqrt{24} = 2 \sqrt{6}$

Thus the normalized ``\underline{w}`` is

$\frac{\underline{w}}{\sqrt{⟨w, w⟩}} = \frac{1}{2 \sqrt{6}} \begin{bmatrix} -2 \\ -2i \\ -4i \end{bmatrix} = \begin{bmatrix} -\frac{1}{\sqrt{6}} \\ -\frac{i}{\sqrt{6}} \\ -\frac{2i}{\sqrt{6}} \end{bmatrix}$
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"

[compat]
QuadGK = "~2.4.2"
StatsPlots = "~0.14.33"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.8.0-beta1"
manifest_format = "2.0"
project_hash = "97dac50cc915a004d479e9b1a154cc224e06420a"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"

[[deps.Arpack]]
deps = ["Arpack_jll", "Libdl", "LinearAlgebra", "Logging"]
git-tree-sha1 = "91ca22c4b8437da89b030f08d71db55a379ce958"
uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
version = "0.5.3"

[[deps.Arpack_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS_jll", "Pkg"]
git-tree-sha1 = "5ba6c757e8feccf03a1554dfaf3e26b3cfc7fd5e"
uuid = "68821587-b530-5797-8361-c406ea357684"
version = "3.5.1+1"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.AxisAlgorithms]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7"
uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
version = "1.0.1"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[deps.Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+1"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.Clustering]]
deps = ["Distances", "LinearAlgebra", "NearestNeighbors", "Printf", "SparseArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "75479b7df4167267d75294d14b58244695beb2ac"
uuid = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
version = "0.14.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.5.0+0"

[[deps.Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.DataValues]]
deps = ["DataValueInterfaces", "Dates"]
git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf"
uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5"
version = "0.4.13"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distances]]
deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.7"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "5a4168170ede913a2cd679e53c2123cb4b889795"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.53"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
version = "1.6.0"

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.8"

[[deps.EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[deps.Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bad72f730e9e91c08d9427d5e8db95478a3c323d"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.4.8+0"

[[deps.FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[deps.FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

[[deps.FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "505876577b5481e50d089c1c68899dfb6faebc62"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.4.6"

[[deps.FFTW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.10+0"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[deps.FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[deps.GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "51d2dfe8e590fbd74e7a842cf6d13d8a2f45dc01"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.6+0"

[[deps.GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "RelocatableFolders", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "df5f5b0450c489fe6ed59a6c0a9804159c22684d"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.64.1"

[[deps.GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "83578392343a7885147726712523c39edc714956"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.64.1+0"

[[deps.GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "83ea630384a13fc4f002b77690bc0afeb4255ac9"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.2"

[[deps.Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[deps.Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "a32d672ac2c967f3deb8a81d828afc739c838a06"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+2"

[[deps.Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[deps.Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[deps.HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[deps.HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+1"

[[deps.HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"]
git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.8"

[[deps.IniFile]]
git-tree-sha1 = "f550e6e32074c939295eb5ea6de31849ac2c9625"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.1"

[[deps.IntelOpenMP_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2018.0.3+2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.Interpolations]]
deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "b15fc0a95c564ca2e0a7ae12c1f095ca848ceb31"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.13.5"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.2+0"

[[deps.KernelDensity]]
deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
git-tree-sha1 = "591e8dc09ad18386189610acafb970032c519707"
uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
version = "0.6.3"

[[deps.LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[deps.LERC_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bf36f528eec6634efc60d7ec062008f171071434"
uuid = "88015f11-f218-50d7-93a8-a6af411a945d"
version = "3.0.0+1"

[[deps.LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[deps.Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "4f00cc36fede3c04b8acf9b2e2763decfdcecfa6"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.13"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
version = "0.6.3"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
version = "7.81.0+0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
version = "1.10.2+0"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+1"

[[deps.Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[deps.Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[deps.Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[deps.Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[deps.Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[deps.Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "c9551dd26e31ab17b86cbd00c2ede019c08758eb"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+1"

[[deps.Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.10"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "e595b205efd49508358f7dc670a940c790204629"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2022.0.0+0"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
version = "2.28.0+0"

[[deps.Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
version = "2022.2.1"

[[deps.MultivariateStats]]
deps = ["Arpack", "LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI", "StatsBase"]
git-tree-sha1 = "7008a3412d823e29d370ddc77411d593bd8a3d03"
uuid = "6f286f6a-111f-5878-ab1e-185364afe411"
version = "0.9.1"

[[deps.NaNMath]]
git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.0"

[[deps.NearestNeighbors]]
deps = ["Distances", "StaticArrays"]
git-tree-sha1 = "ded92de95031d4a8c61dfb6ba9adb6f1d8016ddd"
uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
version = "0.4.10"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
version = "1.2.0"

[[deps.Observables]]
git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
version = "0.4.0"

[[deps.OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.8"

[[deps.Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+1"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
version = "0.3.17+2"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
version = "0.8.1+0"

[[deps.OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ab05aa4cc89736e95915b01e7279e61b1bfe33b8"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.14+0"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
version = "1.8.0"

[[deps.PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[deps.PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "bb16469fd5224100e422f0b027d26c5a25de1200"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.2.0"

[[deps.Plots]]
deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "Unzip"]
git-tree-sha1 = "5f6e1309595e95db24342e56cd4dabd2159e0b79"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.27.3"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Ratios]]
deps = ["Requires"]
git-tree-sha1 = "dc84268fe0e3335a62e315a3a7cf2afa7178a734"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.3"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "dc1e451e15d90347a7decc4221842a022b011714"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.5.2"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.RelocatableFolders]]
deps = ["SHA", "Scratch"]
git-tree-sha1 = "cdbd3b1338c72ce29d9584fdbe9e9b70eeb5adca"
uuid = "05181044-ff0b-4ac5-8273-598c1e38db00"
version = "0.1.3"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
version = "0.7.0"

[[deps.Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "6a2f7d70512d205ca8c7ee31bfa9f142fe74310c"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.12"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "4f6ec5d99a28e1a749559ef7dd518663c5eca3d5"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.3"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "25405d7016a47cf2bd6cd91e66f4de437fd54a07"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.16"

[[deps.StatsPlots]]
deps = ["AbstractFFTs", "Clustering", "DataStructures", "DataValues", "Distributions", "Interpolations", "KernelDensity", "LinearAlgebra", "MultivariateStats", "Observables", "Plots", "RecipesBase", "RecipesPipeline", "Reexport", "StatsBase", "TableOperations", "Tables", "Widgets"]
git-tree-sha1 = "4d9c69d65f1b270ad092de0abe13e859b8c55cad"
uuid = "f3b207a7-027a-5e70-b257-86293d7955fd"
version = "0.14.33"

[[deps.StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "57617b34fa34f91d536eb265df67c2d4519b8b98"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.5"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
version = "1.0.0"

[[deps.TableOperations]]
deps = ["SentinelArrays", "Tables", "Test"]
git-tree-sha1 = "e383c87cf2a1dc41fa30c093b2a19877c83e1bc1"
uuid = "ab02a1b2-a7df-11e8-156e-fb1833f50b87"
version = "1.2.0"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
version = "1.10.0"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.UnicodeFun]]
deps = ["REPL"]
git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
version = "0.4.1"

[[deps.Unzip]]
git-tree-sha1 = "34db80951901073501137bdbc3d5a8e7bbd06670"
uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d"
version = "0.1.2"

[[deps.Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[deps.Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4528479aa01ee1b3b4cd0e6faef0e04cf16466da"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.25.0+0"

[[deps.Widgets]]
deps = ["Colors", "Dates", "Observables", "OrderedCollections"]
git-tree-sha1 = "505c31f585405fc375d99d02588f6ceaba791241"
uuid = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62"
version = "0.6.5"

[[deps.WoodburyMatrices]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3"
uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
version = "0.5.5"

[[deps.XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[deps.XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[deps.Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[deps.Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[deps.Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[deps.Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[deps.Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[deps.Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[deps.Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[deps.Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[deps.Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[deps.Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[deps.Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[deps.Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[deps.Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[deps.Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[deps.Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[deps.Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[deps.Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[deps.Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[deps.Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
version = "1.2.12+1"

[[deps.Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e45044cd873ded54b6a5bac0eb5c971392cf1927"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.2+0"

[[deps.libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
version = "5.0.1+0"

[[deps.libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[deps.libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[deps.libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+1"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
version = "1.41.0+1"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
version = "16.2.1+1"

[[deps.x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[deps.x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"

[[deps.xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╟─3415208a-ab79-11ec-1f49-bdfda5782527
# ╟─e9dba4a1-2b85-4f5d-90fc-87ff6a4317ab
# ╟─bb4bf5a6-c482-42bc-8b4c-e8ae7b48bb1d
# ╟─c9a48dd7-2e94-4e44-80ff-a0b7966f5a2f
# ╟─92098eba-7e9d-4d07-8e38-8599778e4034
# ╟─bef51e5b-b663-4139-8a17-018b473ca188
# ╟─201e9fed-4e7b-480e-a523-e879ab34ec2f
# ╟─50258060-a9b6-468c-a57a-c74c9444fe4a
# ╟─7492da11-24f6-4e04-84a5-831e615c4504
# ╟─36330a8c-e1c6-46ee-8255-b583e5cba6ff
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using PkgSRA, Test, Random
Random.seed!(13)

include("../utils/sindy_turing_utils.jl")
# include("../utils/sparse_regression_functions.jl")
# include("../utils/sindy_statistics_utils.jl")

# Generate test data: Lorenz
include("../examples/example_lorenz.jl")
dat = Array(solve_lorenz_system())
numerical_grad = numerical_derivative(dat, ts)

# Test new interface
alg = slstQuantile(4, 0.2) # Will get extra terms
sindy_library = Dict("cross_terms"=>2,"constant"=>nothing);
test_model = sindyc(dat, numerical_grad, nothing, ts,
                        library=sindy_library,
                        optimizer=alg)

# Built correct terms: names
named_terms = build_term_names(test_model)
@testset "Term Metadata" begin
    @test named_terms ==
        ["x", "y", "z", "", "xx", "xy", "xz", "yy", "yz", "zz"]
end

# Built correct terms: values
library_functions = convert_string2function(test_model.library)
terms = calc_augmented_data(dat, library_functions)

@testset "Data Augmentation" begin
    @test size(terms,1) == length(named_terms)
    @test terms[1:3,:] == dat
    @test terms[4,:] == ones(size(dat,2))
    @test terms[5:10,:] == calc_cross_terms(dat, 2)
end;

# Nonzero term and type tests
t1 = get_nonzero_terms(test_model)
t2 = get_nonzero_terms(core_dyn_true)
# typeof(test_model) <: sindycModel
# @testset "Type Tests" begin
@test typeof(test_model) <: sindycModel
@test issubset(t2, t1)
# end

# Also test cross validation function
null_model = sindyc(dat, 0.1.*randn(size(dat)), nothing, ts,
                        library=sindy_library,
                        optimizer=alg)
# null_model = sindyc(dat, 0.1.*randn(size(dat)), nothing, ts,
#                         library=sindy_library,
#                         use_lasso=true)
err_test = sindy_cross_validate(test_model, dat, numerical_grad)
err_null = sindy_cross_validate(null_model, dat, numerical_grad)

@testset "Cross Validation" begin
    @test err_test < err_null
end

# SINDY tests: up to cubed terms

# Generate test data: Van der Pol
include("../examples/example_vanDerPol.jl")
dat = Array(solve_vdp_system())
numerical_grad = numerical_derivative(dat, ts)

sindy_library = Dict("cross_terms"=>[2,3],
                    "constant"=>nothing);
# test_model2 = sindyc(dat, numerical_grad, nothing, ts,
#                         library=sindy_library,
#                         use_lasso=true,
#                         var_names=["x", "y"])
test_model = sindyc(dat, numerical_grad, nothing, ts,
                        library=sindy_library,
                        optimizer=alg,
                        var_names=["x", "y"])

named_terms = build_term_names(test_model2)
@testset "Cubic term names" begin
    @test named_terms ==
        ["x", "y", "", "xx", "xy", "yy", "xxx", "xxy", "xyy", "yyy"]
end


text:
@testset "Enmap plot" begin
    using Plots; gr()
    using Pixell.Enplot
    register_colors!()
    imap = read_map("data/test.fits", sel=(:,:,1))
    @test (plot(imap, color=:planck); true)
end


text:
using DDF

using Test

#TODO @testset "Meshing"


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

using Markdown
using InteractiveUtils

# ╔═╡ a5df7598-814e-11eb-03c8-b1450e513104
begin
	using Pkg; Pkg.activate("MLJ_env", shared=true)
	using Plots
	using Test
	using LinearAlgebra
	using PlutoUI

	push!(LOAD_PATH, "./src")
	using YaCounter
end

# ╔═╡ 36396ea6-814e-11eb-0d9a-9b4ac175c109
md"""
## Statistics

ref. from book **"Data Science from Scratch"**, Chap 5

$(html"
© Pascal, March 2021
")
"""

# ╔═╡ f1dc5cc6-9018-11eb-1768-7d3c85e942fa
const F = Float64

# ╔═╡ ae9bcffe-8835-11eb-0dd2-c3d39f0ae10e
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ e98ef570-2d8e-457f-af69-e600a879b42e
html"""

"""

# ╔═╡ 0cba28ec-8836-11eb-3dba-df3987947923
html"""

"""

# ╔═╡ c8a24b96-8835-11eb-22f2-954d90b6ca80
md"""
### Describing a Single Set of Data
"""

# ╔═╡ accc7d4e-815b-11eb-03e6-053884a6d2f2
md"""
One obvious description of any dataset is simply the data itself:
"""

# ╔═╡ a5b5ebba-814e-11eb-3945-7d6d1fc6aa1f
begin
	const Num_Friends = Int64[100, 49, 41, 40, 25, 21, 21, 19, 19, 18, 18, 16, 15, 15, 15, 15, 14, 14, 13, 13, 13, 13, 12, 12, 11, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

	const Daily_Minutes = Float64[1., 68.77, 51.25, 52.08, 38.36, 44.54, 57.13, 51.4, 41.42, 31.22, 34.76, 54.01, 38.79, 47.59, 49.1, 27.66, 41.03, 36.73, 48.65, 28.12, 46.62, 35.57, 32.98, 35, 26.07, 23.77, 39.73, 40.57, 31.65, 31.21, 36.32, 20.45, 21.93, 26.02, 27.34, 23.49, 46.94, 30.5, 33.8, 24.23, 21.4, 27.94, 32.24, 40.57, 25.07, 19.42, 22.39, 18.42, 46.96, 23.72, 26.41, 26.97, 36.76, 40.32, 35.02, 29.47, 30.2, 31, 38.11, 38.18, 36.31, 21.03, 30.86, 36.07, 28.66, 29.08, 37.28, 15.28, 24.17, 22.31, 30.17, 25.53, 19.85, 35.37, 44.6, 17.23, 13.47, 26.33, 35.02, 32.09, 24.81, 19.33, 28.77, 24.26, 31.98, 25.73, 24.86, 16.28, 34.51, 15.23, 39.72, 40.8, 26.06, 35.76, 34.76, 16.13, 44.04, 18.03, 19.65, 32.62, 35.59, 39.43, 14.18, 35.24, 40.13, 41.82, 35.45, 36.07, 43.67, 24.61, 20.9, 21.9, 18.79, 27.61, 27.21, 26.61, 29.77, 20.59, 27.53, 13.82, 33.2, 25, 33.1, 36.65, 18.63, 14.87, 22.2, 36.81, 25.53, 24.62, 26.25, 18.21, 28.08, 19.42, 29.79, 32.8, 35.99, 28.32, 27.79, 35.88, 29.06, 36.28, 14.1, 36.63, 37.49, 26.9, 18.58, 38.48, 24.48, 18.95, 33.55, 14.24, 29.04, 32.51, 25.63, 22.22, 19, 32.73, 15.16, 13.9, 27.2, 32.01, 29.27, 33, 13.74, 20.42, 27.32, 18.23, 35.35, 28.48, 9.08, 24.62, 20.12, 35.26, 19.92, 31.02, 16.49, 12.16, 30.7, 31.22, 34.65, 13.13, 27.51, 33.2, 31.57, 14.1, 33.42, 17.44, 10.12, 24.42, 9.82, 23.39, 30.93, 15.03, 21.67, 31.09, 33.29, 22.61, 26.89, 23.48, 8.38, 27.81, 32.35, 23.84]

	const Daily_Hours = Daily_Minutes ./ 60.
end

# ╔═╡ d56f5f5c-816b-11eb-1bf0-bd5824f9bd63
md"""
Let's define some utility functions first
"""

# ╔═╡ e7027f4e-816b-11eb-3919-ffd373f70cc5
sum_of_squares(v::Vector{T}) where T <: Real = sum(v .^ 2)

# ╔═╡ ea1008d4-814e-11eb-3556-8b78bb29f962
friends_count = Counter(Num_Friends; nokey_exception=false, defval=0)

# ╔═╡ 83b0cf5c-814f-11eb-0a3a-9b53f8ea656f
begin
	n = maximum(Num_Friends) + 1
	xs = collect(1:n)
	ys = [friends_count[ix] for ix ∈ xs]

	plot(xs, ys, seriestype = :bar, bins=:scott,
		title="Histogram of Friend Counts",
		lw=1,
		xlims=(0, n),
		ylims=(0, 25),
    	xticks = 0:10:n,
		xlabel="# of friends",
		ylabel="# of people",
		legend=nothing)
end

# ╔═╡ 8265bc76-815b-11eb-1eab-67efd5f15d0e
md"""
Unfortunately, this chart is still too difficult to slip into conversations. So
you start generating some statistics. Probably the simplest statistic is the
number of data points:
"""

# ╔═╡ e89b86b2-815b-11eb-1f9a-7dc9a833233e
num_points = length(Num_Friends)

# ╔═╡ fd2455a0-815b-11eb-222a-7b19df1c0123
md"""
Then largest and smallest values:
"""

# ╔═╡ 095cc486-815c-11eb-0a10-c16a401a3c03
begin
	largest_value, smallest_value = maximum(Num_Friends), minimum(Num_Friends)
end

# ╔═╡ 27f992c2-815c-11eb-31d0-49d06a6b3d39
md"""
Which are just special cases of wanting to know the values in specific positions:
"""

# ╔═╡ 37c684ee-815c-11eb-0710-273f93e7da14
begin
	sorted_values = sort(Num_Friends)

	second_smallest_value = sorted_values[2]
	second_largest_value = sorted_values[end - 1]
	(second_smallest_value, second_largest_value)
end

# ╔═╡ 2d07f7f0-8836-11eb-3fe7-010ea1a62e11
html"""

"""

# ╔═╡ 72f52228-815c-11eb-135e-cd8437424521
md"""
#### Central Tendencies

Usually, we will want some notion of where our data is centered. Most commonly we will use the mean (or average), which is just the sum of the data divided by its count:
"""

# ╔═╡ 8dbdd82a-815c-11eb-3cf9-11b411dc5255
function mean(v::Vector{T})::F where T <: Real
	@assert length(v) > 0
	one(Float64) * sum(v) / length(v)
end

# ╔═╡ c1fbc854-815c-11eb-2e11-37a579e958f1
@test mean(Num_Friends) ≈ 7.333333333333333

# ╔═╡ 16a00e4c-815d-11eb-1785-b965a0359d47
md"""
if you have 10 data points, and you increase the value of any of them by 1, you increase the *mean* by 0.1.
"""

# ╔═╡ 16830248-815d-11eb-118c-c34adc05bdb6
md"""
We will also sometimes be interested in the median, which is the middle-most value (if the number of data points is odd) or the average of the two middle-most values (if the number of data points is even).
"""

# ╔═╡ 166b0b98-815d-11eb-1574-bde5309670cd
function median(v::Vector{T})::T where T <: Real
	@assert length(v) > 0

	n = length(v)
	vₛ = sort(v, rev=false)
 	if isodd(n)
		mix = (1 + n) ÷ 2
		vₛ[mix]
	else
		mix = n ÷ 2
		DT = eltype(vₛ)
		s = vₛ[mix] + vₛ[mix + 1]
		DT <: Integer ? s ÷ DT(2) : s /2.
		# what about Rational and Complex?
	end
end

# ╔═╡ 164e8f4a-815d-11eb-3ef1-7d670678d772
begin
	@test median([1, 10, 2, 9, 5]) == 5
	@test median([1, 9, 2, 10]) == (2 + 9) ÷ 2

	@test median([1., 10., 2., 9., 5.5]) == 5.5
	@test median([1., 9., 2., 10.]) == (2. + 9.) / 2.

	@test median(Num_Friends) == 6
end

# ╔═╡ 15e97c0e-815d-11eb-2207-ab653fc57612
md"""
Clearly, the mean is simpler to compute, and it varies smoothly as our data changes. If we have n data points and one of them increases by some small amount e, then necessarily the mean will increase by e / n. (This makes the mean amenable to all sorts of calculus tricks.)

In order to find the median, however, we have to sort our data (although there are some more efficient technique like binary heap...). And changing one of our data points by a small amount e might increase the median by e, by some number less than e, or not at all (depending on the rest of the data).

However the mean is very sensitive to outliers in our data. If our friendliest user had 200 friends (instead of 100), then the mean would rise to 7.82, while the median would stay the same.
If outliers are likely to be bad data (or otherwise unrepresentative of whatever phenomenon we’re trying to understand), then the mean can sometimes give us a misleading picture.

A generalization of the median is the quantile, which represents the value under which a certain percentile of the data lies (the median represents the value under which 50% of the data lies):
"""

# ╔═╡ 49a1b4a0-8160-11eb-38e3-8dd705fbb346
function quantile(v::Vector{T1}, p::T)::T1 where {T1 <: Real, T <: Real}
	"""Returns the pth-percentile value in x"""
	pix = floor(Int64, p * length(v))
	sort(v, rev=false)[pix]
end

# ╔═╡ 4982ff74-8160-11eb-335d-cd2773020748
begin
	@test quantile(Num_Friends, 0.10) == 1
	@test quantile(Num_Friends, 0.25) == 3
	@test quantile(Num_Friends, 0.5) == median(Num_Friends)
	@test quantile(Num_Friends, 0.75) == 9
	@test quantile(Num_Friends, 0.90) == 13
end

# ╔═╡ 4969769e-8160-11eb-3fa9-4782c4d63ac0
md"""
You might want to look at the mode, or most common value(s):
"""

# ╔═╡ ad124582-8169-11eb-2a3e-f7060063b756
max([1, 2, 2, 3, 3, 4, 5, 6, 6, 7, 6, 3]...)

# ╔═╡ 4935696c-8160-11eb-07ef-059b429dbc9a
function mode(v::Vector{T})::Vector{T} where T <: Real
	"""
	Returns a list, since there might be more than one mode
	"""
	counts = Counter(v)
	max_cnt = values(counts) |> maximum
	[x_i for (x_i, cnt) ∈ counts if cnt == max_cnt]
end

# ╔═╡ d6beae9c-816a-11eb-05fa-092969ea79fd
@test mode(Num_Friends) == Int64[1, 6]

# ╔═╡ 533f5418-8836-11eb-00d4-1d1e2dcb878a
html"""

"""

# ╔═╡ 3c8002a8-816b-11eb-2697-6356a438995d
md"""
#### Dispersion

Dispersion refers to measures of how spread out our data is. Typically they are statistics for which values near zero signify not spread out at all and for which large values (whatever that means) signify very spread out.

For instance, a very simple measure is the range, which is just the difference
between the largest and smallest elements:
"""

# ╔═╡ 3c636418-816b-11eb-2e8d-23bddb884230
range(v::Vector{T}) where T <: Real = maximum(v) - minimum(v)

# ╔═╡ 3c4a8c9a-816b-11eb-24b8-430fadbee06c
@test range(Num_Friends) == 99

# ╔═╡ af74654c-816b-11eb-156c-a343acaed426
md"""
The range is zero precisely when the max and min are equal, which can only happen if the elements of x are all the same, which means the data is not dispersed at all. Conversely, if the range is large, then the max is much larger than the min and the data is more spread out.

Like the median, the range doesn’t really depend on the whole dataset. A dataset whose points are all either 0 or 100 has the same range as a dataset whose values are 0, 100, and lots of 50s. But it seems like the first dataset “should” be more spread out.
"""

# ╔═╡ af5a0c38-816b-11eb-148b-b5964037cba3
function de_mean(v::Vector{T})::Vector{F} where T <: Real
	"""
	Translate v by subtracting its mean (so the result has mean 0)
	"""
	v̄ = mean(v)
	v .- v̄
end

# ╔═╡ af3b67a6-816b-11eb-360c-e5668c8ac456
function variance(v::Vector{T})::F where T <: Real
	"""
	Almost the average squared deviation from the mean
	almost because we use n - 1, to get an unbiased estimate
	"""
	@assert length(v) ≥ 2 "variance requires at least two elements"
	n = length(v)
	dev = de_mean(v)
	one(Float64) * sum_of_squares(dev) / (n - 1)
end

# ╔═╡ 1991ce0a-816d-11eb-3b8f-0fdf9c863bf6
md"""
Now, whatever units our data is in (e.g., “friends”), all of our measures of central tendency are in that same unit. The range will similarly be in that same unit. The variance, on the other hand, has units that are the square of the original units (e.g., “friends squared”). As it can be hard to make sense of these, we often look instead at the *standard deviation*, usually denoted as σ:
"""

# ╔═╡ 28e61e86-816d-11eb-198a-8f9df29bf4d5
σ(v::Vector{T}) where T <: Real = √(variance(v))

# ╔═╡ 28d4fd42-816d-11eb-3ffc-d753aca1fae0
@test 9.02 < σ(Num_Friends) < 9.04

# ╔═╡ 28ba2e5e-816d-11eb-09a4-35e5a7314a2e
md"""
Both the range and the standard deviation have the same outlier problem that we saw earlier for the mean. Using the same example, if our friendliest user had instead 200 friends, the standard deviation would be 14.89 — more than 60% higher!

A more robust alternative computes the difference between the 75th percentile value and the 25th percentile value:
"""

# ╔═╡ 289f7136-816d-11eb-3b5c-e917eef9e556
function interquartile_range(v::Vector{T})::T where T <: Real
	"""
	Returns the difference between the 75%-ile and the 25%-ile
	"""
	quantile(v, 0.75) - quantile(v, 0.25)
end

# ╔═╡ 2885d1fe-816d-11eb-0845-8d7c3aa1e1ba
@test interquartile_range(Num_Friends) == 6

# ╔═╡ ca1377ae-8836-11eb-2366-6355fd427373
html"""

"""

# ╔═╡ 285912ac-816d-11eb-1359-d5ee97622b95
md"""
### Correlation

We’ll first look at covariance, the paired analogue of variance. Whereas variance measures how a single variable deviates from its mean, covariance measures how two variables vary in tandem from their means:
"""

# ╔═╡ 2e768c6a-816e-11eb-303a-f9c19225cadc
function covariance(v₁::Vector{T1}, v₂::Vector{T2})::F where {T1 <: Real, T2 <: Real}
	@assert length(v₁) == length(v₂) "v₁ and v₂ must have same number of elements"
	dot(de_mean(v₁), de_mean(v₂)) / (length(v₁) - 1)
end

# ╔═╡ 2e5a7f2a-816e-11eb-157b-212c9167835d
begin
	@test 22.42 < covariance(Num_Friends, Daily_Minutes) < 22.43
	@test 22.42 / 60. < covariance(Num_Friends, Daily_Hours) < 22.43 / 60
end

# ╔═╡ 2e40e5ba-816e-11eb-38ed-6758053c18da
md"""
Because this number can be hard to interpret, it’s more common to look at the correlation, which divides out the standard deviations of both variables:
"""

# ╔═╡ 2e18ae10-816e-11eb-378e-db14eec0a669
function correlation(v₁::Vector{T1}, v₂::Vector{T2})::F where {T1 <: Real, T2 <: Real}
	"""
	Measures how much v₁ and v₂ vary in tandem about their means
	"""
	σ₁, σ₂ = σ(v₁), σ(v₂)

	σ₁ > zero(F) && σ₂ > zero(F) ? covariance(v₁, v₂) / σ₁ / σ₂ :
		zerot(F)  ## if no variation, correlation is zero
end

# ╔═╡ 2de6a9ce-816e-11eb-2787-6950d6475db8
begin
	@test 0.24 < correlation(Num_Friends, Daily_Minutes) < 0.25
	@test 0.24 < correlation(Num_Friends, Daily_Hours) < 0.25
end

# ╔═╡ 5a29723c-8170-11eb-08a9-71ec55858cc0
plot(Num_Friends, Daily_Minutes, seriestype = :scatter,
	title="Correlation with outlier", lw=1,
	xlims=(0, n), ylims=(0, n), xticks = 0:10:n,
	xlabel="# of friends",
	ylabel="minutes / day",
	legend=nothing)

# ╔═╡ 9c9020e0-8172-11eb-154d-2d7f4d69148c
md"""
The person with 100 friends (who spends only 1 minute per day on the site) is a huge outlier, and correlation can be very sensitive to outliers.
What happens if we ignore him?
"""

# ╔═╡ 5a0aac6a-8170-11eb-188a-bd76519cf2a5
begin
	outlier_ix = findnext(x -> x == 100, Num_Friends, 1) # index of outlier

	num_friends_good = filter(t -> t[1] ≠ outlier_ix,
		enumerate(Num_Friends) |> collect
	) |> a -> map(t -> t[2], a)

	daily_minutes_good =  filter(t -> t[1] ≠ outlier_ix,
		enumerate(Daily_Minutes) |> collect
	) |> a -> map(t -> t[2], a)

	daily_hours_good = daily_minutes_good ./ 60.
end

# ╔═╡ 5014754e-8173-11eb-20ad-83f94f0109ca
begin
	@test 0.57 < correlation(num_friends_good, daily_minutes_good) < 0.58
	@test 0.57 < correlation(num_friends_good, daily_hours_good) < 0.58
end

# ╔═╡ 655d7a4a-8173-11eb-1765-1994d89fcd72
plot(num_friends_good, daily_minutes_good, seriestype = :scatter,
	title="Correlation after removing the outlier", lw=1,
	xlims=(0, n ÷ 2), ylims=(0, n), xticks = 0:10:n,
	xlabel="# of friends",
	ylabel="minutes / day",
	legend=nothing)

# ╔═╡ 8b899ef8-8176-11eb-0d24-3d8972f0e9c4
md"""
Without the outlier, there is a much stronger correlation.
"""

# ╔═╡ Cell order:
# ╟─36396ea6-814e-11eb-0d9a-9b4ac175c109
# ╠═a5df7598-814e-11eb-03c8-b1450e513104
# ╠═f1dc5cc6-9018-11eb-1768-7d3c85e942fa
# ╠═ae9bcffe-8835-11eb-0dd2-c3d39f0ae10e
# ╟─e98ef570-2d8e-457f-af69-e600a879b42e
# ╟─0cba28ec-8836-11eb-3dba-df3987947923
# ╟─c8a24b96-8835-11eb-22f2-954d90b6ca80
# ╟─accc7d4e-815b-11eb-03e6-053884a6d2f2
# ╠═a5b5ebba-814e-11eb-3945-7d6d1fc6aa1f
# ╟─d56f5f5c-816b-11eb-1bf0-bd5824f9bd63
# ╠═e7027f4e-816b-11eb-3919-ffd373f70cc5
# ╠═ea1008d4-814e-11eb-3556-8b78bb29f962
# ╠═83b0cf5c-814f-11eb-0a3a-9b53f8ea656f
# ╟─8265bc76-815b-11eb-1eab-67efd5f15d0e
# ╠═e89b86b2-815b-11eb-1f9a-7dc9a833233e
# ╟─fd2455a0-815b-11eb-222a-7b19df1c0123
# ╠═095cc486-815c-11eb-0a10-c16a401a3c03
# ╟─27f992c2-815c-11eb-31d0-49d06a6b3d39
# ╠═37c684ee-815c-11eb-0710-273f93e7da14
# ╟─2d07f7f0-8836-11eb-3fe7-010ea1a62e11
# ╟─72f52228-815c-11eb-135e-cd8437424521
# ╠═8dbdd82a-815c-11eb-3cf9-11b411dc5255
# ╠═c1fbc854-815c-11eb-2e11-37a579e958f1
# ╟─16a00e4c-815d-11eb-1785-b965a0359d47
# ╟─16830248-815d-11eb-118c-c34adc05bdb6
# ╠═166b0b98-815d-11eb-1574-bde5309670cd
# ╠═164e8f4a-815d-11eb-3ef1-7d670678d772
# ╟─15e97c0e-815d-11eb-2207-ab653fc57612
# ╠═49a1b4a0-8160-11eb-38e3-8dd705fbb346
# ╠═4982ff74-8160-11eb-335d-cd2773020748
# ╟─4969769e-8160-11eb-3fa9-4782c4d63ac0
# ╠═ad124582-8169-11eb-2a3e-f7060063b756
# ╠═4935696c-8160-11eb-07ef-059b429dbc9a
# ╠═d6beae9c-816a-11eb-05fa-092969ea79fd
# ╟─533f5418-8836-11eb-00d4-1d1e2dcb878a
# ╟─3c8002a8-816b-11eb-2697-6356a438995d
# ╠═3c636418-816b-11eb-2e8d-23bddb884230
# ╠═3c4a8c9a-816b-11eb-24b8-430fadbee06c
# ╟─af74654c-816b-11eb-156c-a343acaed426
# ╠═af5a0c38-816b-11eb-148b-b5964037cba3
# ╠═af3b67a6-816b-11eb-360c-e5668c8ac456
# ╟─1991ce0a-816d-11eb-3b8f-0fdf9c863bf6
# ╠═28e61e86-816d-11eb-198a-8f9df29bf4d5
# ╠═28d4fd42-816d-11eb-3ffc-d753aca1fae0
# ╟─28ba2e5e-816d-11eb-09a4-35e5a7314a2e
# ╠═289f7136-816d-11eb-3b5c-e917eef9e556
# ╠═2885d1fe-816d-11eb-0845-8d7c3aa1e1ba
# ╟─ca1377ae-8836-11eb-2366-6355fd427373
# ╟─285912ac-816d-11eb-1359-d5ee97622b95
# ╠═2e768c6a-816e-11eb-303a-f9c19225cadc
# ╠═2e5a7f2a-816e-11eb-157b-212c9167835d
# ╟─2e40e5ba-816e-11eb-38ed-6758053c18da
# ╠═2e18ae10-816e-11eb-378e-db14eec0a669
# ╠═2de6a9ce-816e-11eb-2787-6950d6475db8
# ╠═5a29723c-8170-11eb-08a9-71ec55858cc0
# ╟─9c9020e0-8172-11eb-154d-2d7f4d69148c
# ╠═5a0aac6a-8170-11eb-188a-bd76519cf2a5
# ╠═5014754e-8173-11eb-20ad-83f94f0109ca
# ╠═655d7a4a-8173-11eb-1765-1994d89fcd72
# ╟─8b899ef8-8176-11eb-0d24-3d8972f0e9c4

max_stars_repo_path: Data Science From Scratch/05_Statistics.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
using Documenter, FastLapackInterface

makedocs(sitename = "MyDocumentation")


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_production_vars!(
    model::JuMP.Model,
    g::Unit,
    formulation_prod_vars::Gar1962.ProdVars,
)::Nothing
    prod_above = _init(model, :prod_above)
    segprod = _init(model, :segprod)
    for t in 1:model[:instance].time
        for k in 1:length(g.cost_segments)
            segprod[g.name, t, k] = @variable(model, lower_bound = 0)
        end
        prod_above[g.name, t] = @variable(model, lower_bound = 0)
    end
    return
end

function _add_production_limit_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation_prod_vars::Gar1962.ProdVars,
)::Nothing
    eq_prod_limit = _init(model, :eq_prod_limit)
    is_on = model[:is_on]
    prod_above = model[:prod_above]
    reserve = model[:reserve]
    gn = g.name
    for t in 1:model[:instance].time
        # Objective function terms for production costs
        # Part of (69) of Kneuven et al. (2020) as C^R_g * u_g(t) term
        add_to_expression!(model[:obj], is_on[gn, t], g.min_power_cost[t])

        # Production limit
        # Equation (18) in Kneuven et al. (2020)
        #   as \bar{p}_g(t) \le \bar{P}_g u_g(t)
        # amk: this is a weaker version of (20) and (21) in Kneuven et al. (2020)
        #      but keeping it here in case those are not present
        power_diff = max(g.max_power[t], 0.0) - max(g.min_power[t], 0.0)
        if power_diff < 1e-7
            power_diff = 0.0
        end
        eq_prod_limit[gn, t] = @constraint(
            model,
            prod_above[gn, t] + reserve[gn, t] <= power_diff * is_on[gn, t]
        )
    end
end


text:
module CommunalHelperSyncedZipMoverActivationController

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/SyncedZipMoverActivationController" Controller(
    x::Integer,
    y::Integer,
    colorCode::String="ff00ff",
    zipMoverSpeedMultiplier::Number=1.0,
)

const placements = Ahorn.PlacementDict(
    "Synced Zip Mover Activator (Communal Helper)" => Ahorn.EntityPlacement(
        Controller,
    ),
)

const sprite = "objects/CommunalHelper/syncedZipMoverActivationController/syncedZipMoverActivationController"

function Ahorn.selection(entity::Controller)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Controller)
    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
"""
Describes a river that naturally flows to a reservoir. Its fields correspond
exactly to those of the abstract `River`.
"""
immutable NaturalRiver <: River
  name::AbstractString
  scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1} # 10^6 m^3 / period, indexed with scenarios.
  ext::Dict{AbstractString, Any}

    # Private constructor, without units nor defaut values.
  function NaturalRiver(name::AbstractString,
                        scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1},
                        ext::Dict{AbstractString, Any})
    o = new(name, scenarios, ext)
    _checkConsistency(o)
    return o
  end
end

# Public constructors, with units and keyword arguments.
NaturalRiver(;name::AbstractString=error("Missing reservoir name"),
              scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}=error("Missing scenarios"),
              ext::Dict{AbstractString, Any}=Dict{AbstractString, Any}()) =
  # Always call the internal constructor!
  NaturalRiver(name, scenarios, ext)

copy(r::NaturalRiver; name::AbstractString=r.name, scenarios::Array{TimeSeries.TimeArray{Float64, 1, DateTime, Array{Float64, 1}}, 1}=r.scenarios, ext::Dict{AbstractString, Any}=r.ext) =
  # Always call the internal constructor!
  NaturalRiver(name, scenarios, ext)

Base.isempty(r::NaturalRiver) = false


text:
module Nanoflann

export KDTree, knn

using NearestNeighbors

include("../deps/deps.jl")

type KDTree
    ptr::Ptr{Void}
end

function _free(tree::KDTree)
    ccall((:nanoflann_free_tree, nanoflann_lib), Void, (Ptr{Void},), tree.ptr)
end

function KDTree(data)
    ptr = ccall((:nanoflann_build_tree, nanoflann_lib),
                Ptr{Void}, (Ptr{Cdouble}, Cint, Csize_t),
                data, size(data, 1), size(data, 2))
    tree = KDTree(ptr)
    finalizer(tree, _free)
    tree
end

function NearestNeighbors.knn(tree::KDTree, points, k)
    dim = size(points,1)
    npoints = size(points,2)
    inds   = Array(Csize_t, k, npoints)
    dists2 = Array(Cdouble, k, npoints)
    ccall((:nanoflann_knn, nanoflann_lib),
          Cint, (Ptr{Void}, Ptr{Cdouble}, Csize_t, Csize_t, Cint, Ptr{Csize_t}, Ptr{Cdouble}),
                 tree.ptr,  points,       dim,     npoints, k,    inds,         dists2) == 1 ||
        error("Error in nanoflann_knn")
    (map(Int, inds)+1, dists2)
end

end # module


text:
# Simulation example 1

using CSV
using DataFrames
using Cairo
using Revise
using BayesNegativeBinomial
using BNPRegressionGGA2021
using Distributions
using Gadfly
using Random
using Statistics
using Test
using LinearAlgebra

function simulate_sample(rng, N0, N1, D)
    dϵ = MixtureModel(Normal, [(1.0, 1.0), (-1.0, 0.25)], [0.5, 0.5])
    x0 = [-0.5, 0.5]
    X0 = [ones(N0) rand(rng, x0, N0, D)]
    y0 = X0 * [ones(2); zeros(D - 1)] + rand(rng, dϵ, N0)
    y1 = LinRange(minimum(y0), maximum(y0), N1) |> collect |> x -> repeat(x, 4)
    x1 = kron(x0[[1, 1, 2, 2]], ones(N1))
    x2 = kron(x0[[1, 2, 1, 2]], ones(N1))
    X1 = [ones(4 * N1) x1 x2 zeros(4 * N1, D - 2)]
    return y0, X0, y1, X1
end

# Estimate the density for one realization
p = begin
    N0, N1, D = 1000, 50, 10
    rng = MersenneTwister(1);
    y0, X0, y1, X1 = simulate_sample(rng, N0, N1, D);
    smpl = BNPRegressionGGA2021.NormalSampler(; y0, X0, y1, X1);
    chainf, chainβ = BNPRegressionGGA2021.sample(rng, smpl; mcmcsize = 10000, burnin = 5000);

    fb = mean(chainf);
    plot(x = y0, color = X0[:, 2], Geom.histogram(density = true, bincount = 50))
    p = plot(x = y1, y = fb, color = string.(X1[:, 2]), xgroup = string.(X1[:, 3]),  Geom.subplot_grid(Geom.line))
    img = SVG("figures/fig01.svg", 14cm, 8cm)
    draw(img, p)
    p
end

# Estimate p(\beta[j] | y) for each j, averaged over 100 samples
function get_γchains(R = 1)
    N0, N1, D = 1000, 50, 10
    γchains = zeros(D + 1, 100)
    rng = MersenneTwister(1);
    for i in 1:R
        y0, X0, y1, X1 = simulate_sample(rng, N0, N1, D);
        smpl = BNPRegressionGGA2021.NormalSampler(; y0, X0, y1, X1);
        chainf, chainβ = BNPRegressionGGA2021.sample(rng, smpl; mcmcsize = 5000, burnin = 2500);
        γchains[:, i] .= 1 .- mean([iszero.(chainβ[it]) for it in 1:length(chainβ)])
        println(i)
    end
    return γchains
end
γchains = get_γchains(100)
df = DataFrame(γchains', :auto)
CSV.write("data/data.csv", df)


text:
# Graph analysis tools and support calculation tools

function rangeErrMaxPoint2(fgl1::AbstractDFG, id1, fgl2::AbstractDFG ,id2)
  mv1 = getKDEMax(getBelief(fgl1,id1))
  mv2 = getKDEMax(getBelief(fgl2,id2))
  return norm(mv1[1:2]-mv2[1:2])
end

function rangeCompAllPoses(fgl1::AbstractDFG, fgl2::AbstractDFG)
  ranges = Float64[]
  xx,ll = ls(fgl1)
  for x in xx
    push!(ranges, rangeErrMaxPoint2(fgl1,x,fgl2,x))
  end
  return ranges
end

function rangeCompAllPoses(
    valsbaseline::Dict{Int,Array{Float64,1}},
    fglbaseline::AbstractDFG,
    fgltest::AbstractDFG)

  ranges = Float64[]
  xx,ll = ls(fgltest)
  for x in xx
    mv1 = valsbaseline[fglbaseline.IDs[x]]
    mv2 = getKDEMax(getBelief(fgltest,x))
    push!(ranges, norm(mv1[1:2]-mv2[1:2]))
  end
  return ranges
end

## Compare functions

function compareDensity(a::Normal, b::Normal; tol::Float64=1e-10)
  TP = true
  TP = TP && abs(a.μ - b.μ) < tol
  TP = TP && abs(a.σ - b.σ) < tol
  TP
end

function compareDensity(a::MvNormal, b::MvNormal; tol::Float64=1e-10)::Bool
  TP = true
  TP = TP && norm(a.μ - b.μ)<tol
  TP = TP && sum(norm.(a.Σ.mat - b.Σ.mat))<tol
  return TP
end

#


text:
using JLD
using SharedArrays
using Distributed
include("../examples/solenoid.jl")
function obj_fun(u)
	return u[1]^2 + u[2]^2
end
function obj_fun_erg_avg(s, nSteps)
	nRunup = 500
	u = rand(3)
	J = 0.
	for i = 1:nRunup
		u = next(u,s)
	end
	for i = 1:nSteps
		J += obj_fun(u)/nSteps
		u = next(u,s)
	end
	return J
end
function get_JavgN(sind)
	s = [1.0, 4.0, sind]
	n_rep = 160
	n_pts = 8
	J = zeros(n_pts)
	J_proc = SharedArray{Float64}(n_rep)
	J_proc .= 0.
	Ni = 1200
	N = zeros(Int64, n_pts)
	for i = 1:n_pts
		N[i] = Ni*n_rep
		@show Ni
		t = @distributed for n=1:n_rep
			J_proc[n] = obj_fun_erg_avg(s,Ni)/n_rep
		end
		wait(t)
		J[i] = sum(J_proc)
		Ni = Ni*5
	end
	save("../data/finiteDifference/solenoid/r2_$(sind).jld","s", sind,"J", J, "N", N)
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
import OpenCL; const cl = OpenCL
using Cartesian

macro at_proc(p, ex)
   quote
           remotecall( $p, ()->eval(Main,$(Expr(:quote,ex))))
   end
end

function has64Bit(device)
    amd = "cl_amd_fp64"
    khr = "cl_khr_fp64"
    ext = cl.info(device, :extensions)

    return (khr in ext) || (amd in ext)
end

function countOCL()
    cpu_count = 0
    gpu_count = 0
    for d in cl.devices(:cpu)
        if has64Bit(d)
            cpu_count = 1 # We currently only support one device of each sort.
        end
    end
    for d in cl.devices(:gpu)
        if has64Bit(d)
            gpu_count = 1 # We currently only support one device of each sort.
        end
    end

    return {:gpu => gpu_count, :cpu => cpu_count}
end

function determineWorkload()
    cpus = CPU_CORES
    count = countOCL()

    CPU_OCL = count[:cpu] >= 1
    ncpu_workers =  if CPU_OCL
        iceil((CPU_CORES - count[:gpu] - 1) / 3) # If we have opencl support on the cpu we really don't want to spawn more then one process per 3 cores.
    else
        (CPU_CORES - count[:gpu] - 1) # Native Julia workers
    end

    ngpu_workers = count[:gpu]

    return (ncpu_workers, ngpu_workers, CPU_OCL)
end

function determineCapabilities(cluster :: Bool = false, allow32Bit = false, forceJuliaImpl = false)
    if cluster
        if CPU_OCL
            println("Running simulation on CPU with OpenCL")
            devs = filter(has64Bit, cl.devices(:cpu))
            ctx = cl.Context(devs)
            queue = cl.CmdQueue(ctx)
            return (true, true, ctx, queue)
        elseif GPU_OCL
            println("Running simulation on GPU with OpenCL")
            devs = filter(has64Bit, cl.devices(:gpu))
            ctx = cl.Context(devs)
            queue = cl.CmdQueue(ctx)
            return (true, true, ctx, queue)
        else
            println("Running simulation on CPU")
            return (true, false, nothing, nothing)
        end
    else
        try
            if forceJuliaImpl
                error("forced usage of Julia implementation.")
            end
            if any(map(has64Bit, cl.devices()))
                devs = filter(has64Bit, cl.devices())
                ctx = cl.Context([first(devs)])
                queue = cl.CmdQueue(ctx)
                return (true, true, ctx, queue)
            else
                warn("No OpenCL device with Float64 support found!")
                if allow32Bit
                    warn("Searching for device with Float32 support.")
                    device, ctx, queue = cl.create_compute_context()
                    return (false, true, ctx, queue)
                else
                    throw(Exception())
                end
            end
        catch e
            println("Got exception: $e")
            warn("OpenCL is not supported falling back to Julia computation")
            return (true, false, nothing, nothing)
        end
    end
end

# TODO: find a better way to handle arbitrary dimensions.
function createDisturbance(name :: Symbol, ranges)
    dim = length(ranges)
    if dim == 1
        createDisturbance1(name, ranges)
    elseif dim == 2
        createDisturbance2(name, ranges)
    elseif dim == 3
        createDisturbance3(name, ranges)
    elseif dim == 4
        createDisturbance4(name, ranges)
    else
        error("Can't handle dim $dim")
    end
end

function createDisturbance4(name :: Symbol, ranges)
    vals = Dict[]

    @nloops 4 i d -> ranges[d] begin
        args = @ntuple 4 i
        push!(vals, {100.0 => [name, args...]})
    end

    return vals
end

function createDisturbance3(name :: Symbol, ranges)
    vals = Dict[]

    @nloops 3 i d -> ranges[d] begin
        args = @ntuple 3 i
        push!(vals, {100.0 => [name, args...]})
    end

    return vals
end

function createDisturbance2(name :: Symbol, ranges)
    vals = Dict[]

    @nloops 2 i d -> ranges[d] begin
        args = @ntuple 2 i
        push!(vals, {100.0 => [name, args...]})
    end

    return vals
end

function createDisturbance1(name :: Symbol, range)
    vals = Dict[]

    for i in range
        push!(vals, {100.0 => [name, i]})
    end

    return vals
end

function torange(c :: Dict)
    min = c["min"]
    max = c["max"]
    step = c["step"]
    return min:step:max
end

function torange(c :: Real)
    return c
end

function writedlm_row(io::IO, row, dlm = ',')
    pb = PipeBuffer()
    state = start(row)
    while !done(row, state)
        (x, state) = next(row, state)
        Base.writedlm_cell(pb, x, dlm)
        done(row, state) ? write(pb,'\n') : print(pb,dlm)
    end
    (nb_available(pb) > (16*1024)) && write(io, takebuf_array(pb))
    write(io, takebuf_array(pb))
    nothing
end

function apply_punch_down!(A, x0, y0, a, b)
                d1, d2 = size(A)
                @assert x0 in 1:d1
                @assert y0 in 1:d2
                #dist(x, y) = sqrt((x - x0)^2 + (y-y0)^2)
                #dist(x, y) = abs(x-x0) + abs(y-y0)
                dist(x, y) = floor(sqrt((x - x0)^2 + (y-y0)^2))
                #punch(d) = 1- sech(1/b * d) ^ a
                punch(x)=a*exp(-x^2/(2*b^2))+1
                for j in 1:d2
                    for i in 1:d1
                        d = dist(i, j)
                        A[i,j] *= punch(d)
                    end
                end
            end