text:
"""
Choose a random number from a given distribution as the predictability of a problem.
K: importance of more general information in producing more knowledge. Better be an integer.
max_years: max years a project can take to finish.
"""
function choose_problem(researcher::Researcher, model)
	vurrent_time = model.properties["time"]
	K = model.properties["K"]
	max_years = model.properties["max_years"]
	current_time = model.properties["time"]
	success_probability = researcher.risk_taking
  complexity = 0.0  # not incorporated into the model yet. Can determine time_to_finish. TODO
	time_to_finish = researcher.generality * max_years
	if time_to_finish < 1
		time_to_finish = 1.0
	end
	if K == 0
		generality_importance = 1
	else
		generality_importance = K * time_to_finish
	end
	problem_info = information(success_probability)  * generality_importance # amount of info the project yields
	problem = Problem(success_probability, complexity, time_to_finish, problem_info, current_time)
	return problem
end

"""
    pushproblem!(researcher::Researcher, problem::Problem)

push an already created problem to the history of a researcher
"""
function pushproblem!(researcher::Researcher, problem::Problem)
	push!(researcher.problem_history, problem)
	push!(researcher.publication_citations, 0)
	push!(researcher.publication_success, missing)
end

"
updates researchers profile by one time step (e.g. year)
"
function update_researcher!(researcher::Researcher, model)
	researcher.experience += 1
	if researcher.grant > 0
		researcher.grant -= 1
	end
	finished_problems = Integer[]
	for index in length(researcher.publication_success):-1:1
		if ismissing(researcher.publication_success[index])
			if researcher.problem_history[index].time_to_finish > 1
				researcher.problem_history[index].time_to_finish -= 1
			elseif researcher.problem_history[index].time_to_finish == 1
				researcher.problem_history[index].time_to_finish -= 1
				push!(finished_problems, index)
			end
		end
	end
	for problem_index in finished_problems
		publish!(researcher, problem_index)
	end
	add_citations!(researcher, model)
end

function update_researcher!(model::ABM)
	for researcher in values(model.agents)
		update_researcher!(researcher, model)
	end
end

"""
Adds citations to previous papers
"""
function add_citations!(researcher::Researcher, model::ABM)
	K = model.properties["K"]
	published = findall(a -> !ismissing(a) && a==true, researcher.publication_success)
	for paperind in published
		if K == 0
			generality_importance = 1.0
		else
			generality_importance = K * researcher.received_grants[paperind]
		end
		p1 = 1 - (researcher.problem_history[paperind].success_probability/generality_importance)
		p = p1 < 0 ? 0.0 : p1 > 1 ? 1.0 : p1
		n = round(Int64, 20 -(researcher.problem_history[paperind].start_time - researcher.problem_history[paperind].time_to_finish))  # this makes the number of citations a decreasing function of time until 20 years after which no citations will occur
		if n <= 0
			researcher.publication_citations[paperind] = 0
		else
			d = Binomial(n, p)
			researcher.publication_citations[paperind] += rand(d)
		end
	end
end

"
information(p)

find Shannon's entropy given a probability
"
function information(probability::AbstractFloat)
	info = -log(2, probability)
	return info
end

# TODO: optimize
"""
		information(res::Researcher)

Returns all information produced by a researcher.
"""
function information(res::Researcher)
	return sum([res.problem_history[i].information for i in 1:length(res.publication_success) if .&(!ismissing(res.publication_success[i], res.publication_success[i] == true))])
end

"""
    information(researchers::Array{Researcher})

Returns sum of all information produced by researchers
"""
function information(researchers::Array{Researcher})
	si = 0.0
	for res in researchers
		si += information(res)
	end
	return si
end

"
    publish!(researcher::Researcher, problem_index::Int64)

Check whether a research project is accepted to for publication given its `success_probability`.
"
function publish!(researcher::Researcher, problem_index::Int64)
	p = researcher.problem_history[problem_index].success_probability  # probability to publish, i.e. the project yielded the expected result
	if rand() <= p
		researcher.publication_success[problem_index] = true
	else
		researcher.publication_success[problem_index] = false
	end
end

"""
Years of research grant that researchers can receive.
"""
function get_grant!(model::ABM)
	nresearchers = nagents(model)
	max_years  = model.properties["max_years"]
	ids = collect(allids(model))
	scores = Array{Float64}(undef, nresearchers)
	counter = 1
	for id in ids
		scores[counter] = get_score(model[id], model)
		counter += 1
	end

	# make probabilistic order
	priorities = wsample(ids, scores, nresearchers, replace=false, ordered=false)

	index = 1
	while model.properties["max_grants"] > 0  && index ≤ nresearchers
		resid = priorities[index]
		researcher = model[resid]
		problem = choose_problem(researcher, model)
		asked_years = problem.time_to_finish
		if model.properties["max_grants"] >= asked_years && (asked_years+researcher.grant) <= max_years
			researcher.grant += asked_years
			push!(researcher.received_grants, asked_years)
			model.properties["max_grants"] -= asked_years

			# add the problem to the researchers history
			pushproblem!(researcher, problem)
			# researchers[researcher_index] = researcher
		end
		index += 1
	end
	if model.properties["max_grants"] == 0 && index < nresearchers
		for resid in priorities[index:nresearchers]
			push!(model[resid].received_grants, 0.0)
		end
	end
end

"""
Change a researchers risk taking personality.
"""
function mutate_researcher!(researcher::Researcher, model::ABM)
	μ = model.properties["mutation_rate"]
	risk_range  = model.properties["risk_range"]
	max_years  = model.properties["max_years"]
	risk_mutation = wsample([true, false], [μ, 1-μ])
	years_mutation = wsample([true, false], [μ, 1-μ])
	if risk_mutation
		lbound = researcher.risk_taking-risk_range
		ubound = researcher.risk_taking+risk_range
		if lbound < 0.01
			lbound = 0.01
		end
		if ubound > 0.99
			ubound = 0.99
		end
		researcher.risk_taking = rand(lbound:0.01:ubound)
	end
	if years_mutation
		current = researcher.generality
		lbound = researcher.generality - 1
		ubound = researcher.generality + 1
		lbound < 0.5 ? lbound = 0.5 : lbound
		ubound > max_years ? ubound = max_years : ubound
		researcher.generality = rand(lbound:0.1:ubound)
	end
end

"""
researchers with a min experience can train students.
"""
function train_students!(researcher::Researcher, model::ABM)
  nstudents = 0
  if researcher.experience > model.properties["min_prof_exp"]
    n_students_dist = Poisson(model.properties["average_students"])
    nstudents = rand(n_students_dist)
  else
    return
	end
  if nstudents == 0
    return
  end
  for nstudent in 1:nstudents
    student = Researcher(nextid(model), researcher.risk_taking, researcher.generality)
    mutate_researcher!(student, model)
    add_agent!(student, model)
  end
end

function train_students!(model::ABM)
  for researcher in values(model.agents)
    train_students!(researcher, model)
  end
end

"""
    get_hindex(researcher::Researcher)

Calculates H-index for a researcher.
"""
function get_hindex(researcher::Researcher)
	cited_papers = findall(researcher.publication_citations)
	sort!(cited_papers, rev=true)
	hindex = 0.000001  # I need a positive number for sampling in get_grant!
	for (index, pos) in enumerate(cited_papers)
		if pos >= index
			hindex = index
		else
			break
		end
	end
	return hindex
end

"""
A general function for scoring researchers
"""
function get_score(researcher::Researcher, model::ABM)
	A = model.properties["A"]
	B = model.properties["B"]
	C = model.properties["C"]
	pubs = findall(a -> !ismissing(a) && a==true, researcher.publication_success)
	n_pubs = length(pubs)
	score = 1e-6
	if n_pubs == 0
		return score
	end
	citations = sum(researcher.publication_citations)
	knowledge = sum([researcher.problem_history[i].information for i in pubs])
	exp = researcher.experience
	score += ((n_pubs*A) + (citations*B) + (knowledge*C)) / exp
	return score
end

function retire!(researcher::Researcher, model::ABM)
	if researcher.experience ≥ model.properties["max_exp"]
		kill_agent!(researcher, model)
	end
end

function retire!(model::ABM)
	for researcher in values(model.agents)
		retire!(researcher, model)
	end
end

function needing_grants(model)
	output = Int[]
	for (id, agent) in model.agents
		if agent.grant == 0 && agent.experience >= model.properties["min_prof_exp"]
			push!(output, id)
		end
	end
	return output
end

function having_grants(model)
	output = Int[]
	for (id, agent) in model.agents
		if agent.grant > 0
			push!(output, id)
		end
	end
	return output
end

function researcher_scores(model, need_grants)
	scores = Float64[]
	for resid in need_grants
		score  = get_score(model[resid], model)
		push!(scores, score)
	end
	return scores
end

"""
Only keep the top n researchers each generation, where n=available grant.
Researchers are ranked by the `get_score` function. Excludes those with lowest scores who do not have any ongoing grant.
"""
function exclude_unproductive!(model::ABM)
	# only exclude those with lowest scores who do not have any ongoing grant.
	need_grants = needing_grants(model)
	need_grants_len = length(need_grants)
	popsize = nagents(model)
	available_grants = round(Int, model.properties["max_grants"])

	if need_grants_len == 0 || popsize < available_grants || available_grants > need_grants_len
		return
	end

	# have_grants = having_grants(model)
	scores = researcher_scores(model, need_grants)
	# kill with probability proportional to score
	to_die_count = need_grants_len - available_grants
	to_die = wsample(need_grants, 1 ./ scores, to_die_count, replace=false)
	# sorted_researchers = sortperm(scores, rev=true)
	# to_die = need_grants[sorted_researchers][available_grants+1:end]
	for id in to_die
		kill_agent!(id, model)
	end
end


text:
# This file roughly corresponds to functions documented in the
# Tabs API: https://canvas.instructure.com/doc/api/tabs

"""
    Canvas.tabs(c::Union{Course,Group,User}; kwargs...) -> Vector{Tab}, page_data

Return the paginated list of tabs for the course/group/user.
Return a vector of [`Tab`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List available tabs for a course or group (`GET /api/v1/courses/:course_id/tabs`,
`GET /api/v1/groups/:group_id/tabs`, `GET /api/v1/users/:user_id/tabs`)]
(https://canvas.instructure.com/doc/api/tabs#method.tabs.index)
"""
function tabs(co::Union{Course,Group,User}; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(co))/tabs"; kwargs...)
    return Tab.(json), page_data
end

"""
    Canvas.update_tab(c::Course, t::Tab; kwargs...) -> Tab

Modify an existing tab and return the new details. Return a [`Tab`](@ref).

**Canvas API documentation**:
[Update a tab for a course (`PUT /api/v1/courses/:course_id/tabs/:tab_id`)]
(https://canvas.instructure.com/doc/api/tabs#method.tabs.update)
"""
function update_tab(c::Course, t::Tab; kwargs...)
    json = request("PUT", "/api/v1$(Internals.pid(c))$(Internals.pid(t))"; kwargs...)
    return Tab(json)
end

ext: jl


text:
using Pkg
Pkg.add("Distributions")
Pkg.add("Plots")
using Distributions
using Plots

days = 65

alive_plants = 1000

for i in 1:days
    days = i

    binom = Binomial(1, .999)

    live_dead = rand(binom, alive_plants)

    daily_dead = sum(live_dead[isequal.(live_dead, 0)])

    println("On day $days, $daily_dead plants died")

    alive_plants = alive_plants - daily_dead
end

alive_plants

alive_runs = Vector{Any}()

for j in 1:1000
    days = 65

    alive_plants = 1000

for i in 1:days
    days = i

    binom = Binomial(1, .999)

    live_dead = rand(binom, alive_plants)

    daily_dead = sum(live_dead[isequal.(live_dead, 0)])

    alive_plants = alive_plants - daily_dead
end

push!(alive_runs, alive_plants)

end

histogram(alive_runs)

max_stars_repo_name: saberry/SimOps

max_issues_repo_name: saberry/SimOps

max_forks_repo_name: saberry/SimOps


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

using Markdown
using InteractiveUtils

# ╔═╡ 206337dc-f928-11ea-04ac-278c8fc8079d
using Pkg, DrWatson

# ╔═╡ 20646558-f928-11ea-3351-bf12703651b2
begin
	using StatisticalRethinking
end

# ╔═╡ 0e0d2086-f926-11ea-0470-451b528c9d02
md"## Clip-02-01-02t.jl"

# ╔═╡ 2069098c-f928-11ea-1008-d75cba2186b4
md"## snippet 2.1"

# ╔═╡ 20699ba4-f928-11ea-3f37-8359c913281d
begin
	ways = [0, 3, 8, 9, 0]
	ways = ways / sum(ways)
end

# ╔═╡ 20762c66-f928-11ea-32ec-53742fbaff9b
md"## snippet 2.2"

# ╔═╡ 207be6ea-f928-11ea-32df-cd2ed1f344d1
md"##### With Distributions.jl working with distributions is a little different from R. Instead of having `rbinom`, `dbinom` and `pbinom` we just got the `Binomial` distribution which to work with."

# ╔═╡ aa0f1104-f928-11ea-0830-cbd5ac5bcb63
d = Binomial(9, 0.5)             # Binomial distribution

# ╔═╡ aa0f4d9a-f928-11ea-038e-c31458a61022
d.n                              # Number of trials parameter

# ╔═╡ aa0fd1cc-f928-11ea-1da8-c5f3559d2d0a
rand(d)                          # singe random draw

# ╔═╡ aa1b33a8-f928-11ea-328e-bf25c307896b
rand(d, 10)                      # 10 random draws

# ╔═╡ aa1bc3c0-f928-11ea-2dff-950e37ca54d3
pdf(d, 6)                        # probability density of getting a 6

# ╔═╡ aa261070-f928-11ea-3502-83a01308f7bd
cdf(d, 6)                        # cumulative probability of getting 6

# ╔═╡ aa271b96-f928-11ea-162b-390eddab0c11
pdf(Binomial(9, 0.5), 6)         # probability density of getting a 6

# ╔═╡ 9e7664be-f928-11ea-20b2-653cf2037050
md"## End of clip-02-01-02t.jl"

# ╔═╡ Cell order:
# ╟─0e0d2086-f926-11ea-0470-451b528c9d02
# ╠═206337dc-f928-11ea-04ac-278c8fc8079d
# ╠═20646558-f928-11ea-3351-bf12703651b2
# ╟─2069098c-f928-11ea-1008-d75cba2186b4
# ╠═20699ba4-f928-11ea-3f37-8359c913281d
# ╟─20762c66-f928-11ea-32ec-53742fbaff9b
# ╟─207be6ea-f928-11ea-32df-cd2ed1f344d1
# ╠═aa0f1104-f928-11ea-0830-cbd5ac5bcb63
# ╠═aa0f4d9a-f928-11ea-038e-c31458a61022
# ╠═aa0fd1cc-f928-11ea-1da8-c5f3559d2d0a
# ╠═aa1b33a8-f928-11ea-328e-bf25c307896b
# ╠═aa1bc3c0-f928-11ea-2dff-950e37ca54d3
# ╠═aa261070-f928-11ea-3502-83a01308f7bd
# ╠═aa271b96-f928-11ea-162b-390eddab0c11
# ╟─9e7664be-f928-11ea-20b2-653cf2037050


text:
nums = Set(parse.(Int, readlines("day1/input.txt")))

sums = Dict(i+j => (i, j) for i in nums, j in nums)

for elt in nums
    if (2020-elt) in keys(sums)
        print(elt*prod(sums[2020-elt]))
        break
    end
end


text:
#///////////////////////////////////////
#// File Name: cost_gpu_test.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/01/15
#// Description: Test code for src/cost_gpu.jl
#///////////////////////////////////////

using CUDA
using Distributions
using LinearAlgebra
using Random
using RobotOS

@testset "Cost CUDA Test" begin
prediction_device = "cpu";                                                          # "cpu" or "cuda"
prediction_steps = 12;                                                              # number of steps to look ahead in the future
dto = 0.4;                                                                          # observation update time interval [s]
deterministic = false;
Cep = Matrix(1.0I, 2, 2);                                                           # quadratic position cost matrix
Cu = Matrix(1.0I, 2, 2);                                                            # quadratic control cost matrix
β_pos = 0.6;                                                                        # relative weight between instant and terminal pos cost
β_col = 0.4;                                                                        # relative weight between instant and terminal col cost
α_col = 50.0;                                                                       # scale parameter for exponential collision cost
λ_col = 1.0;                                                                        # bandwidth parameter for exponential collision cost
σ_risk = 0.0;                                                                       # risk-sensitiveness parameteru_norm_max = 2.0;
dtc = 0.02;                                                                         # Euler integration time interval [s]
dtr = 0.1;                                                                          # replanning time interval [s]
dtexec = [0.02, 0.06];                                                              # control insertion duration [s]
tcalc = 0.1;                                                                        # pre-allocated control computation time [s] (< dtr)
u_norm_max = 5.0;                                                                   # maximum control norm [m/s^2]
u_nominal_cand = append!([[0.0, 0.0]],
                         [round.([a*cos(deg2rad(θ)), a*sin(deg2rad(θ))], digits=5)
                          for a = [2., 4.] for θ = 0.:45.:(360. - 45.)])            # nominal control candidate value [ax, ay] [m/s^2]
nominal_search_depth = 1;
# Ego initial state
ego_pos_init_vec = [-5., 0.];                                                       # initial ego position [x, y] [m]
ego_pos_goal_vec = [5., 0.];                                                        # initial ego velocity [vx, vy] [m/s]
# Other parameters
target_speed = 1.0;
sim_horizon = 10.0;                                                                 # simulation time horizon [s]

num_samples = 50;
num_ado_agents = 3;

apos_shared = [0.0, -5.0]
ado_pos_init_dict = Dict{String, Vector{Float64}}();
for ii = 1:num_ado_agents
    ado_pos_init_dict["PEDESTRIAN/$(ii)"] = apos_shared
end
ado_vel_dict = Dict("Any" => MvNormal([0.0, 0.0], Diagonal([0.8, 0.8])));
prediction_rng_seed = 1;

rng = MersenneTwister(prediction_rng_seed);
scene_param = SyntheticSceneParameter(rng);
predictor_param = GaussianPredictorParameter(prediction_steps,
                                             num_samples, deterministic, rng);
cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
cnt_param = ControlParameter(u_norm_max, tcalc, dtexec, dtr,
                             [0.0, 0.0], u_nominal_cand, nominal_search_depth);

scene_loader, controller, w_init, measurement_schedule =
    controller_setup(scene_param,
                     predictor_param,
                     prediction_device=prediction_device,
                     cost_param=cost_param,
                     cnt_param=cnt_param,
                     dtc=dtc,
                     dto=dto,
                     ado_pos_init_dict=ado_pos_init_dict,
                     ado_vel_dict=ado_vel_dict,
                     ego_pos_init_vec=ego_pos_init_vec,
                     ego_vel_init_vec=nothing,
                     ego_pos_goal_vec=ego_pos_goal_vec,
                     target_speed=target_speed,
                     sim_horizon=sim_horizon,
                     verbose=false)

# Nominal Controls
u_array_1 = [[0.0, 0.0] for ii = 1:240];
u_array_2 = [[0.1, 0.05] for ii = 1:240];
u_array = cat(repeat(reshape(u_array_1, (1, length(u_array_1))), inner=(12, 1)),
              repeat(reshape(u_array_2, (1, length(u_array_2))), inner=(12, 1)),
              dims=1);

u_array_gpu_1 = reshape(transpose(hcat(u_array_1...)), 1, 240, 2);
u_array_gpu_1 = cu(repeat(u_array_gpu_1, inner=(12, 1, 1)));
u_array_gpu_2 = reshape(transpose(hcat(u_array_2...)), 1, 240, 2);
u_array_gpu_2 = cu(repeat(u_array_gpu_2, inner=(12, 1, 1)));
u_array_gpu = cat(u_array_gpu_1, u_array_gpu_2, dims=1);

e_state_array_1 = simulate_forward(w_init.e_state, u_array_1, controller.sim_param)
e_state_array_2 = simulate_forward(w_init.e_state, u_array_2, controller.sim_param);
e_state_array = cat(repeat(reshape(e_state_array_1, (1, length(e_state_array_1))), inner=(12, 1)),
                    repeat(reshape(e_state_array_2, (1, length(e_state_array_2))), inner=(12, 1)),
                    dims=1);

ex_array_cpu = Array{Float64, 3}(undef, size(e_state_array, 1), size(e_state_array, 2), 4);
for ii = 1:size(ex_array_cpu, 1)
    for jj = 1:size(ex_array_cpu, 2)
        for kk = 1:size(ex_array_cpu, 3)
            ex_array_cpu[ii, jj, kk] = e_state_array[ii, jj].x[kk]
        end
    end
end
ex_array_gpu = cu(ex_array_cpu);
prediction_horizon = dto*prediction_steps;
target_trajectory = get_nominal_trajectory(w_init.e_state,
                                           ego_pos_goal_vec,
                                           target_speed,
                                           sim_horizon,
                                           prediction_horizon);
target_pos_array = get_target_pos_array(ex_array_gpu, w_init, target_trajectory,
                                        controller.sim_param);
target_pos_array_gpu = cu(target_pos_array);

# Instant Position Cost
inst_pos_cost_array_gpu = instant_position_cost(ex_array_gpu, target_pos_array_gpu,
                                                controller.sim_param.cost_param);
inst_pos_cost_array_cpu = collect(inst_pos_cost_array_gpu);
inst_pos_cost_array =
        map(s -> instant_position_cost(s, target_trajectory, controller.sim_param.cost_param), e_state_array[:, 1:end-1]);
@test all(isapprox.(inst_pos_cost_array_cpu, inst_pos_cost_array, atol=5e-5));

# Instant Control Cost
inst_cnt_cost_array_gpu = instant_control_cost(u_array_gpu, controller.sim_param.cost_param);
inst_cnt_cost_array_cpu = collect(inst_cnt_cost_array_gpu);
inst_cnt_cost_array =
        map(u -> instant_control_cost(u,controller.sim_param.cost_param), u_array);
@test all(isapprox.(inst_cnt_cost_array_cpu, inst_cnt_cost_array, atol=5e-5));

# Instant Collision Cost
prediction_dict = sample_future_ado_positions!(controller.predictor, ado_pos_init_dict);
for key in keys(prediction_dict)
    prediction_dict[key] = repeat(prediction_dict[key], outer=(size(u_array_gpu, 1), 1, 1))
end

ap_array, time_idx_ap_array, control_idx_ex_array, ado_ids_array =
    process_ap_dict(ex_array_gpu, w_init, measurement_schedule, prediction_dict, controller.sim_param);

ap_array_gpu = cu(ap_array);
time_idx_ap_array_gpu = CuArray{Int32, 1}(time_idx_ap_array);
control_idx_ex_array_gpu = CuArray{Int32, 1}(control_idx_ex_array);

inst_col_cost_array_gpu = instant_collision_cost(ex_array_gpu, ap_array_gpu,
                                                 time_idx_ap_array_gpu, control_idx_ex_array_gpu,
                                                 controller.sim_param.cost_param);
inst_col_cost_array_cpu = collect(inst_col_cost_array_gpu);
inst_col_cost_array = similar(inst_col_cost_array_cpu);
for ii = 1:size(inst_col_cost_array, 1)
    for jj = 1:size(inst_col_cost_array, 2)
        for kk = 1:size(inst_col_cost_array, 3)
            ll = control_idx_ex_array[ii];
            mm = time_idx_ap_array[jj];
            inst_col_cost_array[ii, jj, kk] =
                instant_collision_cost(e_state_array[ll, jj], ap_array[ii, mm, kk, :],
                                      controller.sim_param.cost_param);
        end
    end
end
@test all(isapprox.(inst_col_cost_array, inst_col_cost_array_cpu, atol=5e-5))
@test all(isapprox(inst_col_cost_array_cpu[1:50, :, :], inst_col_cost_array_cpu[51:100, :, :], atol=1e-10))
# Instanat Collision Cost with 0 pedestrians
ap_array_gpu_zeros = CuArray{Float32, 4}(undef, size(ap_array_gpu, 1), size(ap_array_gpu, 2), 0, 2);
inst_col_cost_array_gpu_zeros = instant_collision_cost(ex_array_gpu, ap_array_gpu_zeros,
                                                       time_idx_ap_array_gpu, control_idx_ex_array_gpu,
                                                       controller.sim_param.cost_param);
@test size(inst_col_cost_array_gpu_zeros) == (size(inst_col_cost_array_gpu, 1),
                                              size(inst_col_cost_array_gpu, 2),
                                              0);

# Terminal Position Cost
term_pos_cost_array_gpu = terminal_position_cost(ex_array_gpu, target_pos_array_gpu,
                                                 controller.sim_param.cost_param);
term_pos_cost_array_cpu = collect(term_pos_cost_array_gpu);
term_pos_cost_array =
        map(s -> terminal_position_cost(s, target_trajectory, controller.sim_param.cost_param), e_state_array[:, end:end]);
@test all(isapprox.(term_pos_cost_array_cpu, term_pos_cost_array, atol=5e-5))

# Terminal Collision Cost
term_col_cost_array_gpu = terminal_collision_cost(ex_array_gpu, ap_array_gpu,
                                                  time_idx_ap_array_gpu,
                                                  control_idx_ex_array_gpu,
                                                  controller.sim_param.cost_param);
term_col_cost_array_cpu = collect(term_col_cost_array_gpu);
term_col_cost_array = similar(term_col_cost_array_cpu);
for ii = 1:size(inst_col_cost_array, 1)
    for kk = 1:size(inst_col_cost_array, 3)
        ll = control_idx_ex_array[ii];
        mm = time_idx_ap_array[end];
        term_col_cost_array[ii, end, kk] =
            terminal_collision_cost(e_state_array[ll, end], ap_array[ii, mm, kk, :],
                                    controller.sim_param.cost_param);
    end
end
@test all(isapprox.(term_col_cost_array, term_col_cost_array_cpu, atol=5e-5))
@test all(isapprox(term_col_cost_array_cpu[1:50, :, :], term_col_cost_array_cpu[51:100, :, :], atol=1e-10))
# Terminal Collision Cost with 0 pedestrians
term_col_cost_array_gpu_zeros = terminal_collision_cost(ex_array_gpu, ap_array_gpu_zeros,
                                                        time_idx_ap_array_gpu,
                                                        control_idx_ex_array_gpu,
                                                        controller.sim_param.cost_param);
@test size(term_col_cost_array_gpu_zeros) == (size(term_col_cost_array_gpu, 1),
                                              size(term_col_cost_array_gpu, 2),
                                              0);

# Instant Position Cost Gradient
inst_pos_cost_grad_array_gpu = instant_position_cost_gradient(ex_array_gpu[13, :, :],
                                                              target_pos_array_gpu,
                                                              controller.sim_param.cost_param);
inst_pos_cost_grad_array_cpu = collect(inst_pos_cost_grad_array_gpu)

inst_pos_cost_grad_array =
        map(s -> instant_position_cost_gradient(s, target_trajectory, controller.sim_param.cost_param), e_state_array[13, 1:end-1]);
inst_pos_cost_grad_array = Array(transpose(hcat(inst_pos_cost_grad_array...)));
@test all(isapprox.(inst_pos_cost_grad_array_cpu, inst_pos_cost_grad_array, atol=5e-5))

# Instant Collision Cost Gradient
inst_col_cost_grad_array_gpu = instant_collision_cost_gradient(ex_array_gpu[1, :, :],
                                                               ap_array_gpu[1:50, :, :, :],
                                                               time_idx_ap_array_gpu,
                                                               controller.sim_param.cost_param);
inst_col_cost_grad_array_cpu = collect(inst_col_cost_grad_array_gpu);
inst_col_cost_grad_array = similar(inst_col_cost_grad_array_cpu);
for ii = 1:size(inst_col_cost_grad_array, 1)
    for jj = 1:size(inst_col_cost_grad_array, 2)
        for kk = 1:size(inst_col_cost_grad_array, 3)
            mm = time_idx_ap_array[jj];
            inst_col_cost_grad_array[ii, jj, kk, :] =
                instant_collision_cost_gradient(e_state_array[1, jj], ap_array[ii, mm, kk, :],
                                                controller.sim_param.cost_param);
        end
    end
end
@test all(isapprox.(inst_col_cost_grad_array_cpu, inst_col_cost_grad_array, atol=5e-5))
# Instant Collision Cost Gradient with 0 pedestrians
inst_col_cost_grad_array_gpu_zeros = instant_collision_cost_gradient(ex_array_gpu[1, :, :],
                                                                     ap_array_gpu_zeros[1:50, :, :, :],
                                                                     time_idx_ap_array_gpu,
                                                                     controller.sim_param.cost_param);
@test size(inst_col_cost_grad_array_gpu_zeros) == (size(inst_col_cost_grad_array_gpu, 1),
                                                   size(inst_col_cost_grad_array_gpu, 2),
                                                   0,
                                                   size(inst_col_cost_grad_array_gpu, 4));

# Terminal Position Cost Gradient
term_pos_cost_grad_array_gpu = terminal_position_cost_gradient(ex_array_gpu[13, :, :], target_pos_array_gpu,
                                                           controller.sim_param.cost_param);
term_pos_cost_grad_array_cpu = collect(term_pos_cost_grad_array_gpu);
term_pos_cost_grad_array = terminal_position_cost_gradient(e_state_array[13, end], target_trajectory,
                                                           controller.sim_param.cost_param);
term_pos_cost_grad_array = Array(transpose(term_pos_cost_grad_array));
@test all(isapprox(term_pos_cost_grad_array_cpu, term_pos_cost_grad_array, atol=5e-5));

# Terminal Collision Cost Gradient
term_col_cost_grad_array_gpu = terminal_collision_cost_gradient(ex_array_gpu[1, :, :],
                                                                ap_array_gpu[1:50, :, :, :],
                                                                time_idx_ap_array_gpu,
                                                                controller.sim_param.cost_param);
term_col_cost_grad_array_cpu = collect(term_col_cost_grad_array_gpu);
term_col_cost_grad_array = similar(term_col_cost_grad_array_cpu);
for ii = 1:size(inst_col_cost_grad_array, 1)
    for kk = 1:size(inst_col_cost_grad_array, 3)
        mm = time_idx_ap_array[end];
        term_col_cost_grad_array[ii, end, kk, :] =
        terminal_collision_cost_gradient(e_state_array[1, end], ap_array[ii, mm, kk, :],
                                         controller.sim_param.cost_param);
    end
end
@test all(isapprox.(term_col_cost_grad_array_cpu, term_col_cost_grad_array, atol=5e-5))
# Terminal Collision Cost Gradient with 0 pedestrians
term_col_cost_grad_array_gpu_zeros = terminal_collision_cost_gradient(ex_array_gpu[1, :, :],
                                                                      ap_array_gpu_zeros[1:50, :, :, :],
                                                                      time_idx_ap_array_gpu,
                                                                      controller.sim_param.cost_param);
@test size(term_col_cost_grad_array_gpu_zeros) == (size(term_col_cost_grad_array_gpu, 1),
                                                   size(term_col_cost_grad_array_gpu, 2),
                                                   0,
                                                   size(term_col_cost_grad_array_gpu, 4));
end


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

using Markdown
using InteractiveUtils

# ╔═╡ faeb518e-651e-4829-87a9-a2fead7908dd
begin
	N_tr = 100
	x_tr = collect(range(-10.0, 10.0; length=N_tr))
	N_z = 50
	z = collect(range(-10.0, 10.0; length=N_z))
end;

# ╔═╡ 66953fb0-6903-462a-af33-f8db32230cfe
project(x::AbstractVector, κ) = κ * x

# ╔═╡ 11781066-4267-4843-b988-bf90b61574f8
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ 9bd5d65f-4c1b-4248-92df-463432bddafc
begin
	using Pkg
	Pkg.activate(@__DIR__)
	using ApproximateGPs
	using LinearAlgebra
	using Distributions
	using FillArrays
	using Optim
	using ParameterHandling
	using Zygote
	using PDMats
	using Plots
	default(; legend=:outertopright, size=(700, 400))
	ingredients("polyagamma.jl")
	using Random
	rng = Random.MersenneTwister(42);
end

# ╔═╡ ab903524-365d-4529-936f-98ef62b44544
θ_init = (;
    gp = (
        σ² = positive(1.0),
        l = positive(2.0),
    ),
    z,
    jitter = fixed(1e-6),
);

# ╔═╡ 3ccb0a88-63b2-45d8-b723-d613dfe28418
θ_flat_init, unflatten = ParameterHandling.value_flatten(θ_init);

# ╔═╡ e5cd8d0f-9efd-4dff-b97d-63bca5ee4971
function init_var_q(z, ::BernoulliLikelihood, x_tr)
    qω = Vector{PolyaGamma{Int,Float64}}(undef, length(x_tr))
    μ = zeros(length(z))
    Σ = Matrix{Float64}(I(length(z)))
    return μ, Σ, qω
end;

# ╔═╡ 6068ea57-13b0-4a54-a79a-3be978379f76
# Initialize the posterior MvNormal parameters and the posterior distributions of ω
μ, Σ, qω = init_var_q(z, BernoulliLikelihood(), x_tr);

# ╔═╡ 609c8f65-c320-460b-ab81-547f35f74248
build_gp(θ) = GP(θ.σ² * with_lengthscale(SEKernel(), θ.l))

# ╔═╡ 644d6f7d-16a3-4ae6-9463-bc2bacde4485
build_latent_gp(θ) = LatentGP(build_gp(θ.gp), BernoulliLikelihood(), θ.jitter)

# ╔═╡ 1bd89171-b79a-4dfc-bcef-42610f46e366
(f_tr, y_tr) = rand(rng, build_latent_gp(ParameterHandling.value(θ_init))(x_tr));

# ╔═╡ 6e80e323-456f-467a-b1d2-130317f67db2
begin
	p_data = scatter(x_tr, y_tr, label="data", title="Data", xlabel="x")
	plot!(x_tr, f_tr, lw=3.0, label="latent gp")
end

# ╔═╡ 16271c13-8091-4192-895f-3f8afc26fc90
function build_fz_lf(θ)
    lf = build_latent_gp(θ)
    fz = lf(θ.z).fx
    return fz, lf
end

# ╔═╡ f3150669-dbe1-442d-b336-d01cab2babd3
function marginals_to_aug_posterior!(qω::AbstractVector, ps::AbstractVector{<:Normal})
    map!(qω, ps) do p
        PolyaGamma(1, sqrt(abs2(mean(p)) + var(p)))
    end
end

# ╔═╡ 55e30726-d959-4063-8f12-e3cda968f974
project(x::AbstractMatrix, κ) = κ * x * κ'

# ╔═╡ 1b1a79db-af79-4665-a96a-deda8c6633b6
# Update the variational parameters in place
function aug_optimize!(fz::AbstractGPs.FiniteGP, x_tr, y_tr, μ, Σ, qω; niter=10)
    K = ApproximateGPs._chol_cov(fz)
    κ = K \ cov(fz.f, fz.x, x_tr)
    y_sign = sign.(y_tr .- 0.5)
    for _ in 1:niter
        pf = marginals(posterior(SparseVariationalApproximation(Centered(), fz, MvNormal(μ, Σ)))(x_tr))
        marginals_to_aug_posterior!(qω, pf)
        Σ .= Symmetric(inv(inv(K) + project(Diagonal(mean.(qω)), κ)))
        μ .= Σ * (project(y_sign / 2, κ) - K \ mean(fz))
    end
    return μ, Σ, qω
end

# ╔═╡ 5493925c-1cfb-4a87-9468-3379779e56c7
begin
	fz_init, lf_init = build_fz_lf(ParameterHandling.value(θ_init))
	aug_optimize!(fz_init, x_tr, y_tr, μ, Σ, qω)
	u_init = SparseVariationalApproximation(Centered(), fz_init, MvNormal(μ, Σ))
	plot!(p_data, posterior(u_init)(x_tr), label="Initial posterior prediction")
end

# ╔═╡ 2b5e129e-4b80-4b1a-807b-f85004729d48
function kl_term(::BernoulliLikelihood{<:LogisticLink}, aug_variables::AbstractVector{<:PolyaGamma})
    sum(aug_variables) do qω
        c = qω.c
        - abs2(c) * mean(qω) + 2log(cosh(c / 2))
    end
end

# ╔═╡ 71154a8c-6bf2-4b76-b1ca-2f67b90432b6
function _expected_aug_loglik(::BernoulliLikelihood{<:LogisticLink}, qf, qω, y)
    m = mean(qf)
    return  m / 2 - (abs2(m) + var(qf)) * mean(qω) / 2 - log(2)
end

# ╔═╡ 5069065f-2161-42cf-8591-0da4d240f2ed
function expected_aug_loglik(y, q_f, lik, aug_variables)
    return map(y, q_f, aug_variables) do y, q, aug
        _expected_aug_loglik(lik, q, aug, y)
    end
end

# ╔═╡ 15f5b1be-11d0-4b35-bc0f-da9f57d93b6c
function _aug_elbo(sva, fx, y, lik, aug_variables, num_data)
    @assert sva.fz.f === fx.f
    post = posterior(sva)
    q_f = marginals(post(fx.x))
    variational_exp = expected_aug_loglik(y, q_f, lik, aug_variables)

    n_batch = length(y)
    scale = num_data / n_batch
    return sum(variational_exp) * scale - Zygote.@ignore(kl_term(lik, aug_variables)) - ApproximateGPs.kl_term(sva, post)
    # return -ApproximateGPs.kl_term(sva, post)
end

# ╔═╡ 3969259c-9c4f-4497-8f0c-ef658286a535
function aug_elbo(
    sva::SparseVariationalApproximation,
    lfx::AbstractGPs.LatentFiniteGP,
    y::AbstractVector,
    aug_variables;
    num_data=length(y),
)
    @assert sva.fz.f === lfx.fx.f
    return _aug_elbo(sva, lfx.fx, y, lfx.lik, aug_variables, num_data)
end

# ╔═╡ 0b79c288-b640-4ebe-afdf-f85bb0336240
# Compute the loss given theta and the variational parameters while updating mu and Sigma
function loss!(θ, μ, Σ, qω)
    fz, lf = build_fz_lf(θ)
    fx = lf(x_tr)
    Zygote.@ignore aug_optimize!(fz, x_tr, y_tr, μ, Σ, qω)
    augsvgp = SparseVariationalApproximation(Centered(), fz, MvNormal(μ, Σ))
    return -aug_elbo(augsvgp, fx, y_tr, qω)
end

# ╔═╡ 1a2dbd9e-4f92-45b6-b05d-ef845dba79f8
function hp_loss(μ, Σ, qω)
    return θ -> loss!(unflatten(θ), μ, Σ, qω)
end

# ╔═╡ 9024c702-8fb5-43fc-8300-29dd50af699d
hp_loss(μ, Σ, qω)(θ_flat_init)

# ╔═╡ 81df4519-f209-4393-88fc-f369eeab9318
optimisation_result = optimize(
    hp_loss(μ, Σ, qω),
    θ -> only(Zygote.gradient(hp_loss(μ, Σ, qω), θ)),
    θ_flat_init,
    LBFGS(;
        alphaguess=Optim.LineSearches.InitialStatic(; scaled=true),
        linesearch=Optim.LineSearches.BackTracking(),
    ),
    Optim.Options(; iterations=1_000);
    inplace=false,
);

# ╔═╡ b10fd344-a76a-4cc4-957a-adf2befb7144
θ_opt = ParameterHandling.value(unflatten(optimisation_result.minimizer))

# ╔═╡ 8fd6b335-0b90-4534-9785-1ecc8ea698cd
begin
	fz_opt, lf_opt = build_fz_lf(θ_opt)
	aug_optimize!(fz_init, x_tr, y_tr, μ, Σ, qω)
	u_opt = SparseVariationalApproximation(Centered(), fz_opt, MvNormal(μ, Σ))
	plot!(p_data, posterior(u_opt)(x_tr), label="Final posterior prediction")
end

# ╔═╡ Cell order:
# ╠═9bd5d65f-4c1b-4248-92df-463432bddafc
# ╠═faeb518e-651e-4829-87a9-a2fead7908dd
# ╠═ab903524-365d-4529-936f-98ef62b44544
# ╠═1bd89171-b79a-4dfc-bcef-42610f46e366
# ╠═6e80e323-456f-467a-b1d2-130317f67db2
# ╠═6068ea57-13b0-4a54-a79a-3be978379f76
# ╠═5493925c-1cfb-4a87-9468-3379779e56c7
# ╠═3ccb0a88-63b2-45d8-b723-d613dfe28418
# ╠═1a2dbd9e-4f92-45b6-b05d-ef845dba79f8
# ╠═9024c702-8fb5-43fc-8300-29dd50af699d
# ╠═81df4519-f209-4393-88fc-f369eeab9318
# ╠═b10fd344-a76a-4cc4-957a-adf2befb7144
# ╠═8fd6b335-0b90-4534-9785-1ecc8ea698cd
# ╠═e5cd8d0f-9efd-4dff-b97d-63bca5ee4971
# ╠═609c8f65-c320-460b-ab81-547f35f74248
# ╠═644d6f7d-16a3-4ae6-9463-bc2bacde4485
# ╠═16271c13-8091-4192-895f-3f8afc26fc90
# ╠═f3150669-dbe1-442d-b336-d01cab2babd3
# ╠═1b1a79db-af79-4665-a96a-deda8c6633b6
# ╠═66953fb0-6903-462a-af33-f8db32230cfe
# ╠═55e30726-d959-4063-8f12-e3cda968f974
# ╠═0b79c288-b640-4ebe-afdf-f85bb0336240
# ╠═3969259c-9c4f-4497-8f0c-ef658286a535
# ╠═15f5b1be-11d0-4b35-bc0f-da9f57d93b6c
# ╠═2b5e129e-4b80-4b1a-807b-f85004729d48
# ╠═5069065f-2161-42cf-8591-0da4d240f2ed
# ╠═71154a8c-6bf2-4b76-b1ca-2f67b90432b6
# ╟─11781066-4267-4843-b988-bf90b61574f8


text:
const PROPER_SCORING_RULES = "[Gneiting and Raftery (2007), \"Strictly"*
    "Proper Scoring Rules, Prediction, and Estimation\""*
    "](https://doi.org/10.1198/016214506000001437)"
const DOC_FINITE =
    "`AbstractArray{<:Finite}` (multiclass classification)"
const DOC_FINITE_BINARY =
    "`AbstractArray{<:Finite{2}}` (binary classification)"
const DOC_ORDERED_FACTOR =
    "`AbstractArray{<:OrderedFactor}` (classification of ordered target)"
const DOC_ORDERED_FACTOR_BINARY =
    "`AbstractArray{<:OrderedFactor{2}}` "*
    "(binary classification where choice of \"true\" effects the measure)"
const DOC_CONTINUOUS = "`AbstractArray{Continuous}` (regression)"
const DOC_COUNT = "`AbstractArray{Count}`"
const DOC_INFINITE = "AbstractArray{<:Infinite}"
const INVARIANT_LABEL =
    "This metric is invariant to class reordering."
const VARIANT_LABEL =
    "This metric is *not* invariant to class re-ordering"

is_measure_type(::Any) = false

# Each of the following traits, with fallbacks defined in
# StatisticalTraits.jl, make sense for some or all measures:

const MEASURE_TRAITS = [:name,
                        :instances,
                        :human_name,
                        :target_scitype,
                        :supports_weights,
                        :supports_class_weights,
                        :prediction_type,
                        :orientation,
                        :reports_each_observation,
                        :aggregation,
                        :is_feature_dependent,
                        :docstring,
                        :distribution_type]

# # FOR BUILT-IN MEASURES (subtyping Measure)

abstract type Measure <: MLJType end
abstract type Aggregated <: Measure end
abstract type Unaggregated <: Measure end

StatisticalTraits.reports_each_observation(::Type{<:Aggregated}) = false
StatisticalTraits.reports_each_observation(::Type{<:Unaggregated}) = true

# # FALLBACK CHECKS

extra_check(::Measure, args...) = nothing
function _check(measure::Measure, yhat, y)
    check_dimensions(yhat, y)
    extra_check(measure, yhat, y)
end
function _check(measure::Measure, yhat, y, w)
    check_dimensions(yhat, y)
    extra_check(measure, yhat, y, w)
end
function _check(measure::Measure, yhat, y, w::Arr)
    check_dimensions(yhat, y)
    check_dimensions(y, w)
    extra_check(measure, yhat, y, w)
end
function _check(measure::Measure, yhat::Arr{<:UnivariateFinite})
    check_dimensions(yhat, y)
    check_pools(yhat, y)
    extra_check(measure, yhat, y)
end
function _check(measure::Measure,
                yhat::Arr{<:UnivariateFinite},
                y,
                w::Arr)
    check_dimensions(yhat, y)
    check_pools(yhat, y)
    extra_check(measure, yhat, y, w)
end
function _check(measure::Measure,
                yhat::Arr{<:UnivariateFinite},
                y,
                w::AbstractDict)
    check_dimensions(yhat, y)
    check_pools(yhat, y)
    check_pools(yhat, w)
    extra_check(measure, yhat, y, w)
end

# # METHODS TO EVALUATE MEASURES

# See measures/README.md for details

single(::Unaggregated, η̂::Missing, η)          = missing
single(::Unaggregated, η̂,          η::Missing) = missing

const Label = Union{CategoricalValue,Number,AbstractString,Symbol,AbstractChar}

# closure for broadcasting:
single(measure::Measure) = (ηhat, η) -> single(measure, ηhat, η)

call(measure::Unaggregated, yhat, y) = broadcast(single(measure), yhat, y)
function call(measure::Unaggregated, yhat, y, w::Arr)
    unweighted = broadcast(single(measure), yhat, y) # `single` closure below
    return w .* unweighted
end
function call(measure::Unaggregated, yhat, y, weight_given_class::AbstractDict)
    unweighted = broadcast(single(measure), yhat, y) # `single` closure below
    w = @inbounds broadcast(η -> weight_given_class[η], y)
    return w .* unweighted
end

# ## Top level

function (measure::Measure)(args...)
    _check(measure, args...)
    call(measure, args...)
end

# # TRAITS

# user-bespoke measures will subtype `Measure` directly and the
# following will therefore not apply:
StatisticalTraits.supports_weights(::Type{<:Union{Aggregated,Unaggregated}}) =
    true

is_measure_type(::Type{<:Measure}) = true
is_measure(m) = is_measure_type(typeof(m))

# docstring fall-back:
_decorate(s::AbstractString) = "`$s`"
_decorate(v::Vector{<:AbstractString}) = join(_decorate.(v), ", ")
function MMI.docstring(M::Type{<:Measure})
    list = _decorate(instances(M))
    ret = "`$(name(M))` - $(human_name(M)) type"
    isempty(list) || (ret *= " with instances $list")
    ret *= ". "
    return ret
end

# display:
show_as_constructed(::Type{<:Measure}) = true
show_compact(::Type{<:Measure}) = true
Base.show(io::IO, m::Measure) = show(io, MIME("text/plain"), m)

# info (see also src/init.jl):
function ScientificTypes.info(M, ::Val{:measure_type})
    values = Tuple(@eval($trait($M)) for trait in MEASURE_TRAITS)
    return NamedTuple{Tuple(MEASURE_TRAITS)}(values)
end
ScientificTypes.info(m, ::Val{:measure}) = info(typeof(m))

# # AGGREGATION

(::Sum)(v) = sum(skipinvalid(v))
(::Sum)(v::LittleDict) = sum(values(v))

(::Mean)(v) = mean(skipinvalid(v))
(::Mean)(v::LittleDict) = mean(values(v))

(::RootMeanSquare)(v) = sqrt(mean(skipinvalid(v).^2))

aggregate(v, measure) = aggregation(measure)(v)

# aggregation is no-op on scalars:
const MeasureValue = Union{Real,Tuple{<:Real,<:Real}} # number or interval
aggregate(x::MeasureValue, measure) = x

# # UNIVERSAL CALLING SYNTAX

# yhat - predictions (point or probabilisitic)
# X - features
# y - target observations
# w - per-observation weights

function value(measure, yhat, X, y, w)
    vfdep     = Val(is_feature_dependent(measure))
    vsweights = Val(supports_weights(measure) ||
                    supports_class_weights(measure))
    return value(measure, yhat, X, y, w, vfdep, vsweights)
end

# # UNIVERSAL CALLING INTERFACE

#  is feature independent, weights not supported:
value(m, yhat, X, y, w, ::Val{false}, ::Val{false}) = m(yhat, y)

#  is feature dependent:, weights not supported:
value(m, yhat, X, y, w, ::Val{true}, ::Val{false}) = m(yhat, X, y)

#  is feature independent, weights supported:
value(m, yhat, X, y, w,         ::Val{false}, ::Val{true}) = m(yhat, y, w)
value(m, yhat, X, y, ::Nothing, ::Val{false}, ::Val{true}) = m(yhat, y)

#  is feature dependent, weights supported:
value(m, yhat, X, y, w,         ::Val{true}, ::Val{true}) = m(yhat, X, y, w)
value(m, yhat, X, y, ::Nothing, ::Val{true}, ::Val{true}) = m(yhat, X, y)

# # helpers

_scale(x, w::Arr, i) = x*w[i]
_scale(x, ::Nothing, i::Any) = x

function check_pools(ŷ, y)
    levels(y) == levels(ŷ[1]) ||
        error("Conflicting categorical pools found "*
              "in observations and predictions. ")
    return nothing
end

function check_pools(ŷ, w::AbstractDict)
    Set(levels(ŷ[1])) == Set(keys(w)) ||
        error("Conflicting categorical pools found "*
              "in class weights and predictions. ")
    return nothing
end

# # INCLUDE SPECIFIC MEASURES AND TOOLS

include("meta_utilities.jl")
include("roc.jl")
include("confusion_matrix.jl")
include("continuous.jl")
include("finite.jl")
include("probabilistic.jl")
include("loss_functions_interface.jl")

# # DEFAULT MEASURES

default_measure(T, S) = nothing

# Deterministic + Continuous / Count ==> RMS
default_measure(::Type{<:Deterministic},
                ::Type{<:Union{Vec{<:Union{Missing,Continuous}},
                               Vec{<:Union{Missing,Count}}}}) = rms

# Deterministic + Finite ==> Misclassification rate
default_measure(::Type{<:Deterministic},
                ::Type{<:Vec{<:Union{Missing,Finite}}}) = misclassification_rate

# Probabilistic + Finite ==> log loss
default_measure(::Type{<:Probabilistic},
                ::Type{<:Vec{<:Union{Missing,Finite}}}) = log_loss

# Probablistic + Continuous ==> Log loss
default_measure(::Type{<:Probabilistic},
                ::Type{<:Vec{<:Union{Missing,Continuous}}}) = log_loss

# Probablistic + Count ==> Log score
default_measure(::Type{<:Probabilistic},
                ::Type{<:Vec{<:Union{Missing,Count}}}) = log_loss

default_measure(::Type{<:MMI.ProbabilisticDetector},
                ::Type{<:Vec{<:Union{Missing,OrderedFactor{2}}}}) = area_under_curve

default_measure(::Type{<:MMI.DeterministicDetector},
                ::Type{<:Vec{<:Union{Missing,OrderedFactor{2}}}}) = balanced_accuracy

# Fallbacks
default_measure(M::Type{<:Supervised}) = default_measure(M, target_scitype(M))
default_measure(::M) where M <: Supervised = default_measure(M)

default_measure(M::Type{<:Annotator}) = default_measure(M, target_scitype(M))
default_measure(::M) where M <: Annotator = default_measure(M)


text:
using Base.Test
using StrainRecon

@testset "srDistance" begin
(m,n) = (3,2)
M1 = trues(m,n)
M2 = falses(m,n)
w1 = [1.0;0.0]
w2 = [1.0;0.0]
@test abs(srDistance(M1,w1,M2,w2)-sqrt(2*m)) < 1e-12
@test abs(srDistance(M1,w1,M2,w2,p=1) - sqrt(m)) < 1e-12
@test abs(srDistance(M1,w1,M2,w2,p=1,pn=1.0) - m) < 1e-12

M3 = falses(m,n+1)
w3 = [1.0;0.0;0.0]
@test abs(srDistance(M1,w1,M3,w3)-sqrt(2*m)) < 1e-12

w3 = [.5;.5]
M1multi = cat(3,M1,M1)
M2multi = cat(3,M1,M1)
w1multi = cat(2,w1,w1)
w2multi = cat(2,w3,w3)
distances = srDistance(M1multi,w1multi,M2multi,w2multi)
@test abs(distances[2] - sqrt(m/n)) < 1e-12

end


text:
f(v) = let x = 0; for i in v
        while i > 0
            i = max(2, trunc(i / 3)) - 2
            x += i
        end
    end; x
end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
module TestWiringLayers

using Test
using Catlab.WiringDiagrams

# Low-level graph interface
###########################

f = WiringLayer(3, 3)
@test nwires(f) == 0
@test wires(f) == []

add_wire!(f, 1 => 3)
@test has_wire(f, 1 => 3)
@test !has_wire(f, 3 => 1)
@test_throws BoundsError add_wire!(f, 4 => 1) # Source out of bounds
@test_throws BoundsError add_wire!(f, 1 => 4) # Target out of bounds

add_wires!(f, [2 => 2, 3 => 1])
@test nwires(f) == 3
@test wires(f) == [1 => 3, 2 => 2, 3 => 1]

@test nwires(f, 1 => 3) == 1
add_wires!(f, [1 => 3, 1 => 3])
@test nwires(f, 1 => 3) == 3
@test wires(f) == [1 => 3, 1 => 3, 1 => 3, 2 => 2, 3 => 1]
@test out_wires(f, 1) == [1 => 3, 1 => 3, 1 => 3]
rem_wire!(f, 1 => 3)
@test nwires(f, 1 => 3) == 2
rem_wires!(f, 1 => 3)
@test !has_wire(f, 1 => 3)
@test nwires(f, 1 => 3) == 0
@test wires(f) == [2 => 2, 3 => 1]

# High-level categorical interface
##################################

A, B, C = NLayer(1), NLayer(1), NLayer(1)

# Category
@test id(A) == complete_layer(1,1)
@test compose(id(A), id(A)) == id(A)

# Symmetric monoidal category
@test otimes(id(A), id(B)) == id(otimes(A,B))
@test compose(braid(A,B), braid(B,A)) == id(otimes(A,B))

# Diagonals
@test mcopy(A) == complete_layer(1,2)
@test delete(A) == complete_layer(1,0)
@test compose(otimes(mcopy(A), mcopy(B)), otimes(id(A), braid(A,B), id(B))) ==
  mcopy(otimes(A,B))
@test compose(mcopy(A), otimes(mcopy(A), id(A))) == mcopy(A, 3)
@test compose(mcopy(A), otimes(id(A), mcopy(A))) == mcopy(A, 3)
@test compose(mcopy(A), otimes(id(A), delete(A))) == id(A)

# Codiagonals
@test mmerge(A) == complete_layer(2,1)
@test create(A) == complete_layer(0,1)
@test compose(otimes(id(A), braid(B,A), id(B)), otimes(mmerge(A), mmerge(B))) ==
  mmerge(otimes(A,B))
@test compose(otimes(mmerge(A), id(A)), mmerge(A)) == mmerge(A, 3)
@test compose(otimes(id(A), mmerge(A)), mmerge(A)) == mmerge(A, 3)
@test compose(otimes(id(A), create(A)), mmerge(A)) == id(A)

# Wiring diagrams
#################

function roundtrip(layer::WiringLayer)::WiringLayer
  ports = n::Int -> repeat([nothing], n)
  diagram = to_wiring_diagram(layer, ports(layer.ninputs), ports(layer.noutputs))
  wiring_layer_between(diagram, input_id(diagram), output_id(diagram))
end

layer = WiringLayer(3, 3)
add_wires!(layer, [1 => 1, 1 => 2, 2 => 3, 3 => 1, 3 => 3])
@test roundtrip(layer) == layer

end


text:
glsl_version = 130

function __init__()
    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        global glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GLFW.GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        global glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end
end

error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"

function init_renderer(width, height, title::AbstractString)
    # setup GLFW error callback
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(width, height, title)
    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    CImGui.StyleColorsDark()
    # CImGui.StyleColorsClassic()
    # CImGui.StyleColorsLight()

    # setup Platform/Renderer bindings
    CImGui.ImGui_ImplGlfw_InitForOpenGL(window, true)
    CImGui.ImGui_ImplOpenGL3_Init(glsl_version)

    return window, ctx
end

gui_task = nothing

function renderloop(window, ctx, ui=()->nothing, hotloading=false)
    try
        while !GLFW.WindowShouldClose(window)
            GLFW.PollEvents()
            CImGui.ImGui_ImplOpenGL3_NewFrame()
            CImGui.ImGui_ImplGlfw_NewFrame()
            CImGui.NewFrame()

            hotloading ? Base.invokelatest(ui) : ui()

            CImGui.Render()
            GLFW.MakeContextCurrent(window)
            display_w, display_h = GLFW.GetFramebufferSize(window)
            glViewport(0, 0, display_w, display_h)
            glClearColor(0.2, 0.2, 0.2, 1)
            glClear(GL_COLOR_BUFFER_BIT)
            CImGui.ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

            GLFW.MakeContextCurrent(window)
            GLFW.SwapBuffers(window)
            yield()
        end
    catch e
        @error "Error in renderloop!" exception=e
        Base.show_backtrace(stderr, catch_backtrace())
    finally
        CImGui.ImGui_ImplOpenGL3_Shutdown()
        CImGui.ImGui_ImplGlfw_Shutdown()
        CImGui.DestroyContext(ctx)
        GLFW.DestroyWindow(window)
    end
end

function render(ui; width=1280, height=720, title::AbstractString="Demo", hotloading=false)
    window, ctx = init_renderer(width, height, title)
    GC.@preserve window ctx begin
        t = @async renderloop(window, ctx, ui, hotloading)
    end
    return t
end

function gui()
    show_playground_window = true
    show_experiment_window = true
    render_playground_window = playground_window()
    render_experiment_window = experiment_window()
    if gui_task !== nothing
        # ex = InterruptException()
        # Base.throwto(gui_task, ex)
    end
    Plots.theme(:dark)
    global gui_task = Threads.@spawn Threads.@sync render(width = 600, height = 800, title = "MOEA") do
        CImGui.Begin("Start")
        if CImGui.Button("Playground")
            show_playground_window = !show_playground_window
        end
        if show_playground_window
            render_playground_window()
        end
        CImGui.SameLine()
        if CImGui.Button("Experiment")
            show_experiment_window = !show_experiment_window
        end
        if show_experiment_window
            render_experiment_window()
        end
        CImGui.End()
    end
end


text:
##############################################################################
# This files constains utilities to manipulate Julia Expr
##############################################################################

function process_expression!(nonlin_expr::Expr, lin_expr::Dict{Int32, Float64},
                             varlinearities::Vector{Symbol})
    # Get list of all variables in the expression
    extract_variables!(lin_expr, nonlin_expr)
    # Extract linear and constant terms from non-linear expression
    tree = LinearityExpr(nonlin_expr)
    tree = pull_up_constants(tree)
    _, tree, constant = prune_linear_terms!(tree, lin_expr)
    # Make sure all terms remaining in the tree are .nl-compatible
    nonlin_expr = convert_formula(tree)

    # Track which variables appear nonlinearly
    nonlin_vars = Dict{Int32, Float64}()
    extract_variables!(nonlin_vars, nonlin_expr)
    for j in keys(nonlin_vars)
        varlinearities[j] = :Nonlin
    end

    # Remove variables at coeff 0 that aren't also in the nonlinear tree
    for (j, coeff) in lin_expr
        if coeff == 0 && !(j in keys(nonlin_vars))
            delete!(lin_expr, j)
        end
    end

    # Mark constraint as nonlinear if anything is left in the tree
    linearity = nonlin_expr != 0 ? :Nonlin : :Lin

    return nonlin_expr, constant, linearity
end
function process_expression!(nonlin_expr::Real, lin_expr, varlinearities)
    # Special case where body of constraint is constant
    # Return empty nonlinear and linear parts, and use the body as the constant
    0, nonlin_expr, :Lin
end

# For MOI
function process_expression!(nonlin_expr::Expr, lin_expr::Dict{Int32, Float64})
    # Get list of all variables in the expression
    extract_variables!(lin_expr, nonlin_expr)
    # Extract linear and constant terms from non-linear expression
    tree = LinearityExpr(nonlin_expr)
    tree = pull_up_constants(tree)
    _, tree, constant = prune_linear_terms!(tree, lin_expr)
    # Make sure all terms remaining in the tree are .nl-compatible
    nonlin_expr = convert_formula(tree)

    # Track which variables appear nonlinearly
    nonlin_vars = Dict{Int32, Float64}()
    extract_variables!(nonlin_vars, nonlin_expr)

    # Remove variables at coeff 0 that aren't also in the nonlinear tree
    for (j, coeff) in lin_expr
        if abs(coeff) < eps(0.) && !(j in keys(nonlin_vars))
            delete!(lin_expr, j)
        end
    end

    # Mark constraint as nonlinear if anything is left in the tree
    linearity = nonlin_expr != 0 ? :Nonlin : :Lin

    return nonlin_expr, constant, linearity
end

function process_expression!(nonlin_expr::Real, lin_expr)
    # Special case where body of constraint is constant
    # Return empty nonlinear and linear parts, and use the body as the constant
    0, nonlin_expr, :Lin
end

# We need to track linear coeffs of all variables present in the expression tree
extract_variables!(lin_constr::Dict{Int32, Float64}, c) = c
extract_variables!(lin_constr::Dict{Int32, Float64}, c::LinearityExpr) =
    extract_variables!(lin_constr, c.c)
function extract_variables!(lin_constr::Dict{Int32, Float64}, c::Expr)
    if c.head == :ref
        CONFIG[:debug] && println("DEBUG: extract_variables :: variable case: $c")
        if c.args[1] == :x
            varidx = getvidx(c.args[2])
            lin_constr[varidx] = 0
        else
            error("Unrecognized reference expression $c")
        end
    else
        map(arg -> extract_variables!(lin_constr, arg), c.args)
    end
end

add_constant(c, constant::Real) = c + constant
add_constant(c::Expr, constant::Real) = Expr(:call, :+, c, constant)

substitute_vars!(c, x::Array{Float64}) = c
function substitute_vars!(c::Expr, x::Array{Float64})
    if c.head == :ref
        if c.args[1] == :x
            varidx = getvidx(c.args[2])
            c = x[varidx]
        else
            error("Unrecognized reference expression $c")
        end
    else
        if c.head == :call
            # Convert .nl unary minus (:neg) back to :-
            if c.args[1] == :neg
                c.args[1] = :-
            # Convert .nl :sum back to :+
            elseif c.args[1] == :sum
                c.args[1] = :+
            end
        end
        map!(arg -> substitute_vars!(arg, x), c.args, c.args)
    end
    c
end

function evaluate_linear(linear_coeffs::Dict{Int32, Float64}, x::Array{Float64})
    total = 0.0
    for (i, coeff) in linear_coeffs
        total += coeff * x[i]
    end
    total
end

function evaluate_quad(rowidx, colidx, qvals, x::Array{Float64})
   n = length(x)
   mat = sparse(rowidx, colidx, qvals, n, n)
   # This is soooo ugly --xhub
   Q = (mat + mat') - Diagonal(diag(mat))
   total = .5*x'*Q*x
   return total
end


text:
int main() {
        int x;
        int x;
	return 0 ;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
#=
Functionality regarding spatial coordinates of an array, i.e. Longitude and Lattitude
as well as equal area grid handling and anything else related to space.
=#
using Statistics, StatsBase, StaticArrays
export SVector # for equal area grid
#########################################################################
# Spatial indexing
#########################################################################
export spatialidxs, lonlatfirst

"""
    spatialidxs(A::ClimArray) → idxs
Return an iterable that can be used to access all spatial points of `A` with the syntax
```julia
idxs = spatialidxs(A)
for i in idxs
    slice_at_give_space_point = A[i...]
end
```
Works for all types of space (`...` is necessary because `i` is a `Tuple`).
"""
spatialidxs(A::AbDimArray) = spatialidxs(spacestructure(A), A)
function spatialidxs(::LonLatGrid, A)
    lons = (Lon(i) for i in 1:size(A, Lon))
    lats = (Lat(i) for i in 1:size(A, Lat))
    return Iterators.product(lons, lats)
end

"""
    lonlatfirst(A::ClimArray, args...) → B
Permute the dimensions of `A` to make a new array `B` that has first dimension longitude,
second dimension latitude, with the remaining dimensions of `A` following
(useful for most plotting functions). Optional extra dimensions
can be given as `args...`, specifying a specific order for the remaining dimensions.

Example:
```julia
B = lonlatfirst(A)
C = lonlatfirst(A, Time)
```
"""
function lonlatfirst(C, args...)
    permutedims(C, (Lon, Lat, args..., otherdims(C, (Lon, Lat, args...))...))
end

#########################################################################
# Periodicity of longitude
#########################################################################
export lon_distance, wrap_lon, longitude_circshift
"""
    wrap_lon(x)
Wrap given longitude to -180 to 180 degrees.
"""
wrap_lon(x) = @. -180 + (360 + ((x+180) % 360)) % 360

"""
    lon_distance(λ1, λ2, Δλ = 360) → δ
Calculate distance `δ` (also in degrees) between longitudes `λ1, λ2`, but taking into
account the periodic nature of longitude, which has period `Δλ = 360ᵒ`.
"""
function lon_distance(x, y, p = eltype(x)(360))
    moddis = mod(abs(x - y), p)
    min(moddis, p - moddis)
end

"""
    longitude_circshift(X::ClimArray, l = size(X, Lon)÷2, wrap = false)
Perform the same action as `Base.circshift`, but only for the longitudinal dimension
of `X` with shift `l`. If `wrap = true` the longitudes are wrapped to (-180, 180) degrees.
"""
function longitude_circshift(X::ClimArray, l::Int = size(X, Lon)÷2, wrap = false)
    shifts = map(d -> d isa Lon ? l : 0, dims(X))
    shifted_data = circshift(X.data, shifts)
    shifted_lon = mod.(circshift(dims(X, Lon).val, l), 360)
    if wrap; shifted_lon = wrap_lon.(shifted_lon); end
    shifted_lon = vector2range(shifted_lon)
    shifted_dim = Lon(shifted_lon; metadata = dims(X, Lon).metadata)
    new_dims = map(d -> d isa Lon ? shifted_dim : d, dims(X))
    return ClimArray(shifted_data, new_dims; name = X.name, attrib = X.attrib)
end

#########################################################################
# averaging functions over space or time
#########################################################################
using StatsBase

export latmean, spacemean, zonalmean, spaceagg, uniquelats

"""
    zonalmean(A::ClimArray [, W])
Return the zonal mean of `A`. Works for both [`LonLatGrid`](@ref) as well as
[`UnstructuredGrid`](@ref). Optionally provide statistical weights `W`.
These can be the same `size` as `A` or only having the same latitude structure as `A`.
"""
zonalmean(A::AbDimArray, W = nothing) = zonalmean(spacestructure(A), A, W)
zonalmean(::LonLatGrid, A::AbDimArray, W) = dropagg(mean, A, Lon, W)

"""
    latmean(A::ClimArray)
Return the latitude-mean `A` (mean across dimension `Lat`).
This function properly weights by the cosine of the latitude.
"""
function latmean(A::AbDimArray)
    lw = _latweights(A)
    if ndims(A) > 1
        return dropagg(sum, dimwise(*, A, lw), Lat)
    else
        return sum(A .* lw)
    end
end

# Warning!!! `_latweights` divides by the weight sum, because it is intended to be
# used only with the `sum` function (for A)
_latweights(A::AbDimArray) = _latweights(dims(A, Lat))
function _latweights(A::Lat)
    we = cosd.(A.val)
    we ./= sum(we)
    return ClimArray(we, (A,))
end

using StatsBase

"""
    spacemean(A::ClimArray [, W]) = spaceagg(mean, A, W)
Average given `A` over its spatial coordinates.
Optionally provide statistical weights in `W`.
"""
spacemean(A, exw=nothing) = spaceagg(mean, A, exw)

"""
    spaceagg(f, A::ClimArray, W = nothing)
Aggregate `A` using function `f` (e.g. `mean, std`) over all available space (i.e.
longitude and latitude) of `A`, weighting every part of `A` by its spatial area.

`W` can be extra weights, to weight each spatial point with. `W` can either be a
`ClimArray` with same spatial information as `A`, or having exactly same dimensions as `A`.
"""
spaceagg(f, A::AbDimArray, exw=nothing) = spaceagg(spacestructure(A), f, A, exw)
function spaceagg(::LonLatGrid, f, A::AbDimArray, w=nothing)
    wtype = spaceweightassert(A, w)
    cosweights = repeat(cosd.(dims(A, Lat).val)', size(A, Lon))
    if dimindex(A, Lon) > dimindex(A, Lat)
        error("At the moment this function assumes that Lon preceeds Lat, use `permutdims`.")
    end
    other = otherdims(A, (Lon, Lat))
    # pre-calculate weights if possible
    if wtype == :no
        W = weights(cosweights)
    elseif wtype == :d2
        W = weights(cosweights .* Array(w))
    end
    # if A is only 2-dimensional, operation is trivial and single number
    if ndims(A) == 2
        return f(A, W)
    end
    # do the weighted average
    oidxs = otheridxs(A, (Lon(), Lat()))
    if wtype != :dany
        r = map(i -> f(view(A, i), W), oidxs)
    else
        # TODO: This multiplication .* here assumes that the lon-lat grid is the
        # first dimension.
        r = map(i -> f(view(A, i), weights(view(w, i) .* cosweights)), oidxs)
    end
    return ClimArray(r, other, A.name)
end

function spaceweightassert(A, w)
    if !isnothing(w)
        wdims = dims(w)
        # Check that longitude/latitude match even if the actual dimensions are not
        # identical (because they are e.g. not loaded from the same file)
        if length(wdims) == 2
            @assert Lat ∈ basetypeof.(wdims)
            @assert Lon ∈ basetypeof.(wdims)
            @assert val.(wdims) == val.(dims(A, (Lon, Lat)))
            wtype = :d2
        else
            @assert basetypeof.(wdims) == basetypeof.(dims(A))
            @assert size(w) == size(A)
            wtype = :dany
        end
    else
        wtype = :no
    end
    return wtype
end

#########################################################################
# Hemispheric sum/difference
#########################################################################
export hemispheric_means, hemispheric_functions

"""
    hemispheric_functions(A::ClimArray) → north, south
Return two arrays `north, south`, by splitting `A` to its northern and southern hemispheres,
appropriately translating the latitudes of `south` so that both arrays have the same
latitudinal dimension (and thus can be compared and do opperations between them).
"""
hemispheric_functions(A) = hemispheric_functions(spacestructure(A), A)
function hemispheric_functions(::LonLatGrid, A)
    nh = A[Lat(Between(0,  90))]
    sh = A[Lat(Between(-90, 0))]
    # TODO: this can be a function "reverse dim"
    di = dimindex(sh, Lat)
    newdims = [dims(sh)...]
    newdims[di] = dims(nh, Lat)
    data = reverse(Array(sh); dims = di)
    sh = ClimArray(data, (newdims...,))
    return nh, sh
end

"""
    hemispheric_means(A) → nh, sh
Return the (proper) averages of `A` over the northern and southern hemispheres.
Notice that this function explicitly does both zonal as well as meridional averaging.
Use [`hemispheric_functions`](@ref) to just split `A` into two hemispheres.
"""
hemispheric_means(A, args...) = hemispheric_means(spacestructure(A), A, args...)
function hemispheric_means(::LonLatGrid, A::AbDimArray)
    @assert hasdim(A, Lat)
    if hasdim(A, Lon)
        B = zonalmean(A)
    else
        B = A
    end
    nh = latmean(B[Lat(Between(0,  90))])
    sh = latmean(B[Lat(Between(-90, 0))])
    return nh, sh
end

latitudes(A) = latitudes(spacestructure(A), A)
latitudes(::LonLatGrid, A) = dims(A, Lat).val


text:
######################################################################
# df.jl
#
# Discrete filter
######################################################################

struct DF{OL, S<:Sensor} <: AbstractFilter
    b::Matrix{Float64}      # the actual discrete belief
    n::Int64                # number of cells per side
    cell_size::Float64      # width of each cell, in meters
    sensor::S               # sensor model used in filtering
    obs_list::OL            # list of observations
end
function DF(m::SearchDomain, n::Int64, sensor::Sensor, obs_list=0:0)
    b = ones(n, n) / (n * n)
    return DF(b, n, m.length/n, sensor, obs_list)
end

function update!(df::DF, p::Pose, o)
    n = df.n
    bp_sum = 0.0

    for txi = 1:n, tyi = 1:n
        if df.b[txi, tyi] > 0.0
            tx = (txi-0.5) * df.cell_size
            ty = (tyi-0.5) * df.cell_size
            df.b[txi, tyi] *= O(df.sensor, (tx,ty,0.0,0.0), p, o)
            bp_sum += df.b[txi, tyi]
        end
    end

    # normalize
    for txi = 1:n, tyi = 1:n
        df.b[txi, tyi] /= bp_sum
    end
end

# returns x, y value
# Assumes that the belief sums to one
function centroid(df::DF)
    x_val = 0.0; y_val = 0.0
    for x = 1:df.n, y = 1:df.n
        x_val += (x-.5) * df.b[x,y]
        y_val += (y-.5) * df.b[x,y]
    end
    return x_val*df.cell_size, y_val*df.cell_size
end

function centroid(d::Matrix{Float64}, L::Float64)
    x_val = 0.0; y_val = 0.0
    n = size(d,1)
    cell_size = L / n
    d_sum = 0.0
    for x = 1:n, y = 1:n
        x_val += (x-.5) * d[x,y]
        y_val += (y-.5) * d[x,y]
        d_sum += d[x,y]
    end
    return x_val*cell_size/d_sum, y_val*cell_size/d_sum
end

function covariance(df::DF)
    mu_x = mu_y = 0.0
    c_xx = c_xy = c_yy = 0.0
    for xi = 1:df.n, yi = 1:df.n
        x = (xi-0.5)*df.cell_size
        y = (yi-0.5)*df.cell_size

        mu_x += x * df.b[xi,yi]
        mu_y += y * df.b[xi,yi]

        c_xx += df.b[xi,yi] * x * x
        c_yy += df.b[xi,yi] * y * y
        c_xy += df.b[xi,yi] * x * y
    end
    c_xx -= (mu_x * mu_x)
    c_yy -= (mu_y * mu_y)
    c_xy -= (mu_x * mu_y)
    return [c_xx+1e-4 c_xy; c_xy c_yy+1e-4]
end

# add some to ensure we are
function covariance(d::Matrix{Float64}, L::Float64)
	mu_x, mu_y = centroid(d, L)
	c_xx = c_xy = c_yy = 0.0
	n = size(d,1)
	cell_size = L / n
	d_sum = 0.0
	for xi = 1:n
		for yi = 1:n
			x = (xi-0.5)*cell_size
			y = (yi-0.5)*cell_size

			c_xx += d[xi,yi] * x * x
			c_yy += d[xi,yi] * y * y
			c_xy += d[xi,yi] * (x - mu_x) * (y - mu_y)
			d_sum += d[xi,yi]
		end
	end
	c_xx = c_xx/d_sum - (mu_x * mu_x)
	c_yy = c_yy/d_sum - (mu_y * mu_y)
	c_xy = c_xy/d_sum

	# add 1e-3 to ensure we are positive definite
	return [c_xx+1e-3 c_xy; c_xy c_yy+1e-3]
end

# Returns the entropy of the distribution.
entropy(df::DF) = entropy(df.b)

# I had this before... a normalized version?
#entropy(df::DF) = entropy(df.b) / log(df.n * df.n)

reset!(f::DF) = fill!(f.b, 1.0/(f.n*f.n))

function print_belief(df::DF)
	for y = df.n:-1:1
		for x = 1:(df.n-1)
			@printf "%.2f," df.b[x,y]
		end
		@printf "%.2f\n" df.b[df.n,y]
	end
	#@printf "%.2f",
end

# I never use these...
function noiseless(x::Vehicle, theta::LocTuple)
	noiseless(x, x.sensor, theta)
end
function noiseless(x::Vehicle, s::BearingOnly, theta::LocTuple)
	true_bearing(x, theta)
end
function noiseless(x::Vehicle, s::DirOmni, theta::LocTuple)
	rel_bearing = x.heading - true_bearing(x, theta)
	if rel_bearing < 0.0
		rel_bearing += 360.0
	end
	rel_int = round(Int, rel_bearing, RoundDown) + 1

	return s.means[rel_int]
end


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, Base.Test
prob = deepcopy(prob_ode_2Dlinear)
prob.tspan = (1.0,0.0)

sol = solve(prob,DP5(),dt=-1/4)

if !is_windows()
  sol2 = solve(prob,dopri5(),dt=-1/4)
  @test sol.t ≈ sol2.t
end

sol = solve(prob,DP5(),dt=-1/4,tstops=[0.5])
@test sol.t == [1.0,.75,.5,0]

sol = solve(prob,RK4(),dt=-1/4)
@test sol.t == [1.00,.75,.5,.25,0]

sol = solve(prob,RK4(),dt=-1/4,tstops=[0.5,0.33])
@test sol.t ≈ [1.00,.75,.5,.33,0.08,0]


text:
module Patterns
    using ..Util

    include("pattern.jl")
    include("toexpr.jl")
    include("patterns_syntax.jl")
    include("exports.jl")
end


text:
matrix = Array{Float64,2}(100,100)
matrix[31,42] = pi

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TravelerProblem

include("Test.jl");
include("Tools.jl");
include("BruteForce.jl")
include("LinearProgramming.jl")
include("CheapestRutes.jl")

end # module


text:
# This file includes code that was formerly a part of Julia. License is MIT: http://julialang.org/license

UTF32String(data::Vector{Char}) = UTF32String(reinterpret(UInt32, data))

# UTF-32 basic functions
next(s::UTF32String, i::Int) = (Char(s.data[i]), i+1)
lastindex(s::UTF32String) = length(s.data) - 1
length(s::UTF32String) = length(s.data) - 1

codeunit(s::UTF32String) = UInt32
ncodeunits(s::UTF32String) = length(s.data)

if isdefined(Base, :iterate)
    function iterate(s::UTF32String, i::Int = firstindex(s))
        i > length(s) && return nothing
        return next(s, i)
    end
end

reverse(s::UTF32String) = UTF32String(reverse!(copy(s.data), 1, length(s)))

sizeof(s::UTF32String) = sizeof(s.data) - sizeof(UInt32)

const empty_utf32 = UTF32String(UInt32[0])

convert(::Type{UTF32String}, c::Char) = UTF32String(UInt32[c, 0])
convert(::Type{UTF32String}, s::UTF32String) = s

function convert(::Type{UTF32String}, str::AbstractString)
    len, flags = unsafe_checkstring(str)
    buf = Vector{UInt32}(undef, len+1)
    out = 0
    @inbounds for ch in str ; buf[out += 1] = ch ; end
    @inbounds buf[out + 1] = 0 # NULL termination
    UTF32String(buf)
end

function convert(::Type{UTF8String},  str::UTF32String)
    dat = str.data
    len = sizeof(dat) >>> 2
    # handle zero length string quickly
    len <= 1 && return empty_utf8
    # get number of bytes to allocate
    len, flags, num4byte, num3byte, num2byte = unsafe_checkstring(dat, 1, len-1)
    flags == 0 && @inbounds return UTF8String(copyto!(Vector{UInt8}(undef, len), 1, dat, 1, len))
    return encode_to_utf8(UInt32, dat, len + num2byte + num3byte*2 + num4byte*3)
end

function convert(::Type{UTF32String}, str::UTF8String)
    dat = str.data
    # handle zero length string quickly
    sizeof(dat) == 0 && return empty_utf32
    # Validate UTF-8 encoding, and get number of words to create
    len, flags = unsafe_checkstring(dat)
    # Optimize case where no characters > 0x7f
    flags == 0 && @inbounds return fast_utf_copy(UTF32String, UInt32, len, dat, true)
    # has multi-byte UTF-8 sequences
    buf = Vector{UInt32}(undef, len+1)
    @inbounds buf[len+1] = 0 # NULL termination
    local ch::UInt32, surr::UInt32
    out = 0
    pos = 0
    @inbounds while out < len
        ch = dat[pos += 1]
        # Handle ASCII characters
        if ch <= 0x7f
            buf[out += 1] = ch
        # Handle range 0x80-0x7ff
        elseif ch < 0xe0
            buf[out += 1] = ((ch & 0x1f) << 6) | (dat[pos += 1] & 0x3f)
        # Handle range 0x800-0xffff
        elseif ch < 0xf0
            pos += 2
            ch = get_utf8_3byte(dat, pos, ch)
            # Handle surrogate pairs (should have been encoded in 4 bytes)
            if is_surrogate_lead(ch)
                # Build up 32-bit character from ch and trailing surrogate in next 3 bytes
                pos += 3
                surr = ((UInt32(dat[pos-2] & 0xf) << 12)
                        | (UInt32(dat[pos-1] & 0x3f) << 6)
                        | (dat[pos] & 0x3f))
                ch = get_supplementary(ch, surr)
            end
            buf[out += 1] = ch
        # Handle range 0x10000-0x10ffff
        else
            pos += 3
            buf[out += 1] = get_utf8_4byte(dat, pos, ch)
        end
    end
    UTF32String(buf)
end

function convert(::Type{UTF32String}, str::UTF16String)
    dat = str.data
    len = sizeof(dat)
    # handle zero length string quickly (account for trailing \0)
    len <= 2 && return empty_utf32
    # get number of words to create
    len, flags, num4byte = unsafe_checkstring(dat, 1, len>>>1)
    # No surrogate pairs, do optimized copy
    (flags & UTF_UNICODE4) == 0 && @inbounds return UTF32String(copyto!(Vector{UInt32}(undef, len), dat))
    local ch::UInt32
    buf = Vector{UInt32}(undef, len)
    out = 0
    pos = 0
    @inbounds while out < len
        ch = dat[pos += 1]
        # check for surrogate pair
        if is_surrogate_lead(ch) ; ch = get_supplementary(ch, dat[pos += 1]) ; end
        buf[out += 1] = ch
    end
    UTF32String(buf)
end

function convert(::Type{UTF16String}, str::UTF32String)
    dat = str.data
    len = sizeof(dat)
    # handle zero length string quickly
    len <= 4 && return empty_utf16
    # get number of words to allocate
    len, flags, num4byte = unsafe_checkstring(dat, 1, len>>>2)
    # optimized path, no surrogates
    num4byte == 0 && @inbounds return UTF16String(copyto!(Vector{UInt16}(undef, len), dat))
    return encode_to_utf16(dat, len + num4byte)
end

function convert(::Type{UTF32String}, str::ASCIIString)
    dat = str.data
    @inbounds return fast_utf_copy(UTF32String, UInt32, length(dat), dat, true)
end

function convert(::Type{UTF32String}, dat::AbstractVector{UInt32})
    @inbounds return fast_utf_copy(UTF32String, UInt32, length(dat), dat, true)
end

convert(::Type{UTF32String}, data::AbstractVector{Int32}) =
    convert(UTF32String, reinterpret(UInt32, convert(Vector{Int32}, data)))

convert(::Type{UTF32String}, data::AbstractVector{Char}) =
    convert(UTF32String, map(UInt32, data))

convert(::Type{T}, v::AbstractVector{S}) where {T<:AbstractString, S<:Union{UInt32,Char,Int32}} =
    convert(T, utf32(v))

# specialize for performance reasons:
function convert(::Type{T}, data::AbstractVector{S}) where {T<:ByteString, S<:Union{UInt32,Char,Int32}}
    s = IOBuffer(Vector{UInt8}(undef, length(data)), read=true, write=true)
    truncate(s,0)
    for x in data
        print(s, Char(x))
    end
    convert(T, String(take!(s)))
end

convert(::Type{Vector{UInt32}}, str::UTF32String) = str.data
convert(::Type{Array{UInt32}},  str::UTF32String) = str.data

unsafe_convert(::Type{Ptr{T}}, s::UTF32String) where {T<:Union{UInt32,Int32,Char}} =
    convert(Ptr{T}, pointer(s))

function convert(T::Type{UTF32String}, bytes::AbstractArray{UInt8})
    isempty(bytes) && return empty_utf32
    nb = length(bytes)
    nb & 3 != 0 && throw(UnicodeError(UTF_ERR_ODD_BYTES_32,0,0))
    b1 = bytes[1]
    b2 = bytes[2]
    b3 = bytes[3]
    b4 = bytes[4]
    if b1 == 0 && b2 == 0 && b3 == 0xfe && b4 == 0xff
        offset = 1
        swap = false
    elseif b1 == 0xff && b2 == 0xfe && b3 == 0 && b4 == 0
        offset = 1
        swap = true
    else
        offset = 0
        swap = false
    end
    len = nb ÷ 4 - offset
    d = Vector{UInt32}(undef, len + 1)
    if swap
        @inbounds for i in 1:len
            ib = i + offset
            b1 = UInt32(bytes[ib * 2 - 1])
            b2 = UInt32(bytes[ib * 2])
            b3 = UInt32(bytes[ib * 2 + 1])
            b4 = UInt32(bytes[ib * 2 + 2])
            d[i] = (b1 << 24) | (b2 << 16) | (b3 << 8) | b4
        end
    else
        unsafe_copyto!(Ptr{UInt8}(pointer(d)), pointer(bytes, offset * 4 + 1), len * 4)
    end
    d[end] = 0 # NULL terminate
    UTF32String(d)
end

function isvalid(::Type{UTF32String}, str::Union{Vector{UInt32}, Vector{Char}})
    for c in str
        @inbounds if !isvalid(Char, UInt32(c)) ; return false ; end
    end
    return true
end
isvalid(str::Vector{Char}) = isvalid(UTF32String, str)

utf32(x) = convert(UTF32String, x)

utf32(p::Ptr{UInt32}, len::Integer) = utf32(unsafe_wrap(Array, p, len))
utf32(p::Union{Ptr{Char}, Ptr{Int32}}, len::Integer) = utf32(convert(Ptr{UInt32}, p), len)
function utf32(p::Union{Ptr{UInt32}, Ptr{Char}, Ptr{Int32}})
    len = 0
    while unsafe_load(p, len+1) != 0; len += 1; end
    utf32(p, len)
end

function map(f, s::UTF32String)
    d = s.data
    out = similar(d)
    out[end] = 0

    @inbounds for i = 1:(length(d)-1)
        c2 = f(Char(d[i]))
        if !isa(c2, Char)
            throw(UnicodeError(UTF_ERR_MAP_CHAR, 0, 0))
        end
        out[i] = (c2::Char)
    end
    UTF32String(out)
end

# Definitions for C compatible strings, that don't allow embedded
# '\0', and which are terminated by a '\0'
containsnul(s::ByteString) = containsnul(unsafe_convert(Ptr{Cchar}, s), sizeof(s))
containsnul(s::Union{UTF16String,UTF32String}) = findfirst(isequal(0), s.data) != length(s.data)

if sizeof(Cwchar_t) == 2
    const WString = UTF16String
    const wstring = utf16
elseif sizeof(Cwchar_t) == 4
    const WString = UTF32String
    const wstring = utf32
end
wstring(s::Cwstring) = wstring(convert(Ptr{Cwchar_t}, s))

# Cwstring is defined in c.jl, but conversion needs to be defined here
# to have WString
function unsafe_convert(::Type{Cwstring}, s::WString)
    if containsnul(s)
        throw(ArgumentError("embedded NUL chars are not allowed in C strings: $(repr(s))"))
    end
    return Cwstring(unsafe_convert(Ptr{Cwchar_t}, s))
end

# pointer conversions of ASCII/UTF8/UTF16/UTF32 strings:
pointer(x::Union{ByteString,UTF16String,UTF32String}) = pointer(x.data)
pointer(x::ByteString, i::Integer) = pointer(x.data)+(i-1)
pointer(x::Union{UTF16String,UTF32String}, i::Integer) = pointer(x)+(i-1)*sizeof(eltype(x.data))

# pointer conversions of SubString of ASCII/UTF8/UTF16/UTF32:
pointer(x::SubString{T}) where {T<:ByteString} = pointer(x.string.data) + x.offset
pointer(x::SubString{T}, i::Integer) where {T<:ByteString} = pointer(x.string.data) + x.offset + (i-1)
pointer(x::SubString{T}) where {T<:Union{UTF16String,UTF32String}} = pointer(x.string.data) + x.offset*sizeof(eltype(x.string.data))
pointer(x::SubString{T}, i::Integer) where {T<:Union{UTF16String,UTF32String}} = pointer(x.string.data) + (x.offset + (i-1))*sizeof(eltype(x.string.data))

"""
    utf32(s)

Create a UTF-32 string from a byte array, array of `Char` or `UInt32`, or any other string
type. (Conversions of byte arrays check for a byte-order marker in the first four bytes, and
do not include it in the resulting string.)

Note that the resulting `UTF32String` data is terminated by the NUL codepoint (32-bit zero),
which is not treated as a character in the string (so that it is mostly invisible in Julia);
this allows the string to be passed directly to external functions requiring NUL-terminated
data. This NUL is appended automatically by the `utf32(s)` conversion function. If you have
a `Char` or `UInt32` array `A` that is already NUL-terminated UTF-32 data, then you can
instead use `UTF32String(A)` to construct the string without making a copy of the data and
treating the NUL as a terminator rather than as part of the string.
"""
utf32(s)

"""
    utf32(::Union{Ptr{Char},Ptr{UInt32},Ptr{Int32}} [, length])

Create a string from the address of a NUL-terminated UTF-32 string. A copy is made; the
pointer can be safely freed. If `length` is specified, the string does not have to be
NUL-terminated.
"""
utf32(::Union{Ptr{Char},Ptr{UInt32},Ptr{Int32}}, length=length)


text:
using LibCURL
using Test

# Just testing loading of the library and a simple library call.

@testset "LibCURL" begin
    @testset "test escape" begin
        curl = curl_easy_init()
        curl == C_NULL && error("curl_easy_init() failed")

        function testescape(s, esc_s)
            b_arr = curl_easy_escape(curl, s, sizeof(s))
            us = unsafe_string(b_arr)
            us != esc_s && error("escaping $s failed")
            curl_free(b_arr)
            @test us == esc_s
        end

        testescape("hello world", "hello%20world")
        testescape("hello %world", "hello%20%25world")
        testescape("hello%world", "hello%25world")

        curl_easy_cleanup(curl)
    end

    @testset "download" begin
        curl = curl_easy_init()
        curl == C_NULL && error("curl_easy_init() failed")

        url = "https://github.com/JuliaWeb/LibCURL.jl/blob/master/README.md"
        filename = joinpath(@__DIR__, "README.md")

        @testset "julia download" begin
            try
                download(url, filename)
                @test isfile(filename)
            finally
                rm(filename; force=true)
            end
        end

        @testset "system curl" begin
            try
                run(`curl -g -L -f -o $filename $url`)
                @test isfile(filename)
            finally
                rm(filename; force=true)
            end
        end

        curl_easy_cleanup(curl)
    end

    include("ssl.jl")
end


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 49e8e6f6-0483-11eb-0a74-d94aa536c9ab
using LinearAlgebra, Plots, Zygote, PlutoUI

# ╔═╡ fccca5fc-0597-11eb-1033-a5cc3b71d854
md"""
# Constrained convex optimization

**Michiel Stock**
STMO

## Equality constraints

Consider the following optimization problem:

$$\min_{\mathbf{x}} f(\mathbf{x})$$
$$\text{subject to } g(\mathbf{x})=0\,.$$

For every point $\mathbf{x}$ on the surface $g(\mathbf{x})=0$, the gradient $\nabla g(\mathbf{x})$ is normal to this surface. This can be shown by considering a point $\mathbf{x}+\boldsymbol{\epsilon}$, also on the surface. If we make a Taylor expansion around $\mathbf{x}$, we have

$$g(\mathbf{x}+\boldsymbol{\epsilon})\approx g(\mathbf{x}) + \boldsymbol{\epsilon}^\top\nabla g(\mathbf{x})\,.$$

Given that both $\mathbf{x}$ and $\mathbf{x}+\boldsymbol{\epsilon}$ lie on the surface it follows that $g(\mathbf{x}+\boldsymbol{\epsilon})= g(\mathbf{x})$. In the limit that $||\boldsymbol{\epsilon}||\rightarrow 0$ we have that $\boldsymbol{\epsilon}^\top\nabla g(\mathbf{x})=0$. Because $\boldsymbol{\epsilon}$ is parallel to the surface $g(\mathbf{x})$, it follows that $\nabla g(\mathbf{x})$ is normal to the surface.

We seek a point $\mathbf{x}^\star$ on the surface such that $f(\mathbf{x})$ is minimized. For such a point, it should hold that the gradient w.r.t. $f$ should be parallel to $\nabla g$. Otherwise, it would be possible to give a small 'nudge' to $\mathbf{x}^\star$ in the direction of $\nabla f$ to decrease the function value, which would indicate that $\mathbf{x}^\star$ is not a minimizer. This figures below illustrate this point.

$$\nabla f(\mathbf{x}^\star) + \nu \nabla g (\mathbf{x}^\star)=0\,,$$
with $\nu\neq 0$ called the *Lagrange multiplier*. The constrained minimization problem can also be represented by a *Lagrangian*:
$$L(\mathbf{x}, \nu) 	\equiv f(\mathbf{x}) + \nu g(\mathbf{x})\,.$$
The constrained stationary condition is obtained by setting $\nabla_\mathbf{x} L(\mathbf{x}, \nu) =0$, the condition $\partial  L(\mathbf{x}, \nu)/\partial \nu=0$ leads to the constraint equation $g(\mathbf{x})=0$.

"""

# ╔═╡ 5e4eabba-0483-11eb-09a4-31ae6c582d3e
f((x1, x2)) = 2x1^2 +1.4x1*x2 + x2^2 - 0.3x1 + x1

# ╔═╡ 470048ec-0598-11eb-0d48-67357ed1f12c
md"In addition to a $g(\mathbf{x})$, we also implement a parametric version where $g(\mathbf{x})=0$ for plotting purposes."

# ╔═╡ 5503c842-0591-11eb-0104-359b9ebd5e98
md"Show the countour of $f(\mathbf{x})$:"

# ╔═╡ 47067fc8-058c-11eb-224c-6fe2e470107f
@bind show_f_contour CheckBox()

# ╔═╡ c369c65e-0598-11eb-0a90-dd4db070c99b
show_f_contour

# ╔═╡ 60328136-0591-11eb-1929-b7515a8c03d9
md"Show the countour of $g(\mathbf{x})$:"

# ╔═╡ d10b14c2-058c-11eb-3a43-7bca5d042969
@bind show_g_contour CheckBox()

# ╔═╡ ba835a96-0598-11eb-1d54-9365e3ea337f
show_g_contour

# ╔═╡ 681a031a-0591-11eb-2fec-8d94eca758f8
md"Show the function $g(\mathbf{x})=0$:"

# ╔═╡ 13467b42-058d-11eb-3023-b9d3f145335b
@bind show_g0_constraint CheckBox()

# ╔═╡ ca227ac2-0598-11eb-3a39-f1a90769ac52
show_g0_constraint

# ╔═╡ d8276f26-0598-11eb-1626-95f87cf6de22
md"Show the gradients:"

# ╔═╡ 60c41848-0592-11eb-39ef-7103ae834990
@bind show_gradients CheckBox()

# ╔═╡ d16723ca-0598-11eb-02f4-f1c5d45095a5
show_gradients

# ╔═╡ 95dc760a-058e-11eb-2149-898b2776079f
@bind t Slider(0:0.01:2π, default=4.3)

# ╔═╡ 41d8c8d0-0591-11eb-2cd7-3b53f91ef589
t

# ╔═╡ aef30ee8-0597-11eb-2203-ed0f9b14f0e4
t

# ╔═╡ 00a009a2-0599-11eb-341f-3364198e1753
md"We plot the objective value on the equality constraint."

# ╔═╡ 13b4ad4a-0599-11eb-0043-6bfc0eb23389
md"At the minimizer, the gradients are parallel."

# ╔═╡ 7e77e628-0707-11eb-1b2f-27b5c5819260
# compute angle between two vectors
angle(u, v) = acos(dot(u, v) / (norm(u) * norm(v)))

# ╔═╡ 358f21ca-0599-11eb-3bfb-0bb12a1465bb
md"""
## Inequality constraints

The same argument can be made for inequality constraints, i.e. solving

$$\min_{\mathbf{x}} f(\mathbf{x})$$
$$\text{subject to } g(\mathbf{x})\leq0\,.$$

Here, two situations can arise:

- **Inactive constraint**: the minimizer of $f$ lies in the region where $g(\mathbf{x}) < 0$. This corresponds to a Lagrange multiplier $\nu=0$. Note that the solution would be the same if the constraint was not present.
- **Active constraint**: the minimizer of $f$ lies in the region where $g(\mathbf{x}) > 0$. The solution of the constrained problem will lie on the bound where $g(\mathbf{x})=0$, similar to the equality-constrained problem and corresponds to a Lagrange multiplier $\nu>0$.

For both cases, the product $\nu g(\mathbf{x})=0$, the solution should thus satisfy the following conditions:
$$g(\mathbf{x}) \leq 0$$
$$\nu \geq 0$$
$$\nu g(\mathbf{x})=0\,.$$
These are called the *Karush-Kuhn-Tucker* conditions.

It is relatively straightforward to extend this framework towards multiple constraints (equality and inequality) by using several Lagrange multipliers.
"""

# ╔═╡ 2d7ae208-0599-11eb-1a91-277b764ad45e
md"Change the location of the constraint"

# ╔═╡ c98b4f9a-0592-11eb-2034-5f05e792cbd2
@bind R Slider(10:0.5:25, default=16)

# ╔═╡ bd7b55ae-0599-11eb-1056-ed3682fe8160
g((x1, x2)) = (x1-10)^2 + (x2-17)^2 - R^2

# ╔═╡ c10838a2-0599-11eb-32e3-1579a926f9b7
g0(t) = [R*cos(t) + 10, R*sin(t) + 17]

# ╔═╡ d2c4ffe2-0485-11eb-0e20-dd52446fa3a6
x = g0(t)

# ╔═╡ 6e89822c-070b-11eb-06f3-8b8e20aaa57d
f(x)

# ╔═╡ e944567a-0596-11eb-2cf1-75f196bab4d1
f'(x) / norm(f'(x))

# ╔═╡ 02673cec-0597-11eb-2db4-e9f7d9047ac0
g'(x) / norm(g'(x))

# ╔═╡ c04a0ffe-0707-11eb-2a74-b3b258b4aec6
angle(f'(x), g'(x))  # close to 0 or pi if stationary point

# ╔═╡ f47ce686-0485-11eb-2a33-4bd0b150b18b
g(x)

# ╔═╡ b0cf9a40-058f-11eb-13dc-b7c35328ef2a
begin
	plot(t -> f(g0(t)), 0:0.01:2π, label="f(x) where g(x)=0")
	xlabel!("t")
	ylabel!("objective")
	scatter!([t], [f(g0(t))], color=:purple, label="x")
end

# ╔═╡ be35a06c-0596-11eb-2988-61a79c7e6a87
R

# ╔═╡ a08e7230-0591-11eb-3956-fb6c2d756494
fconstr(x) = g(x) ≤ 0.0 ? f(x) : NaN

# ╔═╡ 81629936-0593-11eb-3fca-afa68e2ace3d
@bind t_approx Slider(0.1:0.05:10.0, default=1.0)

# ╔═╡ a1402cbc-0593-11eb-33b4-09a0c654eefd
t_approx

# ╔═╡ cc4209ee-0593-11eb-3762-e76539e2a455
Î₋(u) = u < 0.0 ? -(1/t_approx)* log(-u) : Inf

# ╔═╡ a398740e-0593-11eb-2f81-c92861e080b8
begin
	plot(-3:0.001:2, Î₋, label="logarithmic barier (t=$t_approx)")
	plot!([-3, 0, 0], [0, 0, 5], ls=:dash, label="indicator function")
	ylims!(-2, 5)
	xlabel!("u")
end

# ╔═╡ dce25666-0594-11eb-24f9-85eba4949d1e
fsoft(x) = g(x) ≤ 0.0 ? t_approx * f(x) + Î₋(g(x)) : NaN

# ╔═╡ 701efa2a-0596-11eb-2592-39eced2fb875
@bind show_soft_grads CheckBox()

# ╔═╡ 7c5d381a-0596-11eb-1dfe-0b7074093cd6
show_soft_grads

# ╔═╡ 7c934678-0599-11eb-1e6b-1f2c9003509a
md"## Utilities"

# ╔═╡ 865bbb8e-0483-11eb-0622-b9a8a6ba5d4c
x1min, x1max = -8.0, 15.0

# ╔═╡ 9736f388-0483-11eb-069a-0d5ed87efb73
x2min, x2max = -8.0, 12.0

# ╔═╡ f39f0106-0483-11eb-3617-5fb04388489e
grads(f, x) = 0.1f'(x) |> df -> Tuple([e] for e in df)

# ╔═╡ 828cbe72-0483-11eb-33b2-7f4f37db889f
begin
	plot(colorbar=false)
	show_f_contour && contourf!(x1min:0.1:x1max, x2min:0.1:x2max, (x1, x2) -> f((x1, x2)), label="f(x)", colorbar=true, color=:Blues, aspect=:equal)
	show_g_contour && contour!(x1min:0.1:x1max, x2min:0.1:x2max, (x1, x2) -> g((x1, x2)), color=:Reds)
	show_g0_constraint && plot!([g0(t)[1] for t in 0:0.01:2π], [g0(t)[2] for t in 0:0.01:2π], lw=2, color=:red, label="g(x)=0")
	show_g0_constraint && scatter!([x[1]], [x[2]], label="x", color=:purple)
	scatter!([0], [0], label="x* (no constraint)", color=:green)
	show_gradients && quiver!([x[1]], [x[2]], quiver=grads(f, x), color=:blue)
	show_gradients && quiver!([x[1]], [x[2]], quiver=grads(g, x), color=:orange)
	xlims!(x1min, x1max)
	ylims!(x2min, x2max)
	xlabel!("x1")
	ylabel!("x2")
end

# ╔═╡ c2640938-0591-11eb-058c-5949d1e316d3
begin
	contourf(x1min:0.01:x1max, x2min:0.01:x2max, (x1, x2) -> fconstr((x1, x2)), 			color=:Blues)
	show_g_contour && contour!(x1min:0.1:x1max, x2min:0.1:x2max, (x1, x2) -> g((x1, x2)), color=:Reds)
	contour!(x1min:0.1:x1max, x2min:0.1:x2max, (x1, x2) -> f((x1, x2)), 			color=:Blues)
	scatter!([x[1]], [x[2]], label="x", color=:purple)
	scatter!([0], [0], label="x* (no constraint)", color=:green)
	show_g0_constraint && plot!([g0(t)[1] for t in 0:0.01:2π], [g0(t)[2] for t in 0:0.01:2π], lw=2, color=:red, label="g(x)=0")
	show_gradients && quiver!([x[1]], [x[2]], quiver=grads(f, x), color=:blue)
	show_gradients && quiver!([x[1]], [x[2]], quiver=grads(g, x), color=:orange)
	xlims!(x1min, x1max)
	ylims!(x2min, x2max)
end

# ╔═╡ 52bacf14-0595-11eb-33c9-dd5c40125360
∇fsoft(x) = g(x) ≤ 0.0 ? fsoft'(x) : [0.0, 0.0]

# ╔═╡ 0b93638a-0595-11eb-3096-af3cc7202aa5
begin
	contourf(x1min:0.01:x1max, x2min:0.01:x2max, (x1, x2) -> fconstr((x1, x2)), 			color=:Blues)

	contour!(x1min:0.1:x1max, x2min:0.1:x2max, (x1, x2) -> f((x1, x2)), 			color=:Blues)
	show_g_contour && contour!(x1min:0.1:x1max, x2min:0.1:x2max, (x1, x2) -> fsoft((x1, x2)), color=:Reds)
	scatter!([x[1]], [x[2]], label="x", color=:purple)
	scatter!([0], [0], label="x* (no constraint)", color=:green)
	show_g0_constraint && plot!([g0(t)[1] for t in 0:0.01:2π], [g0(t)[2] for t in 0:0.01:2π], lw=2, color=:red, label="g(x)=0")
	show_soft_grads && quiver!(x1min:2:x1max, (x2min:2:x2max)', quiver=(x1, x2)->-0.05.*∇fsoft((x1, x2)), color=:pink)
	xlims!(x1min, x1max)
	ylims!(x2min, x2max)
end

# ╔═╡ 0a5783e2-0596-11eb-1e38-a5da40325159

# ╔═╡ Cell order:
# ╠═49e8e6f6-0483-11eb-0a74-d94aa536c9ab
# ╟─fccca5fc-0597-11eb-1033-a5cc3b71d854
# ╠═5e4eabba-0483-11eb-09a4-31ae6c582d3e
# ╠═bd7b55ae-0599-11eb-1056-ed3682fe8160
# ╟─470048ec-0598-11eb-0d48-67357ed1f12c
# ╠═c10838a2-0599-11eb-32e3-1579a926f9b7
# ╠═41d8c8d0-0591-11eb-2cd7-3b53f91ef589
# ╠═d2c4ffe2-0485-11eb-0e20-dd52446fa3a6
# ╠═6e89822c-070b-11eb-06f3-8b8e20aaa57d
# ╠═f47ce686-0485-11eb-2a33-4bd0b150b18b
# ╟─5503c842-0591-11eb-0104-359b9ebd5e98
# ╟─47067fc8-058c-11eb-224c-6fe2e470107f
# ╟─c369c65e-0598-11eb-0a90-dd4db070c99b
# ╟─60328136-0591-11eb-1929-b7515a8c03d9
# ╟─d10b14c2-058c-11eb-3a43-7bca5d042969
# ╟─ba835a96-0598-11eb-1d54-9365e3ea337f
# ╟─681a031a-0591-11eb-2fec-8d94eca758f8
# ╟─13467b42-058d-11eb-3023-b9d3f145335b
# ╟─ca227ac2-0598-11eb-3a39-f1a90769ac52
# ╟─d8276f26-0598-11eb-1626-95f87cf6de22
# ╟─60c41848-0592-11eb-39ef-7103ae834990
# ╟─d16723ca-0598-11eb-02f4-f1c5d45095a5
# ╟─828cbe72-0483-11eb-33b2-7f4f37db889f
# ╠═95dc760a-058e-11eb-2149-898b2776079f
# ╠═aef30ee8-0597-11eb-2203-ed0f9b14f0e4
# ╟─00a009a2-0599-11eb-341f-3364198e1753
# ╟─b0cf9a40-058f-11eb-13dc-b7c35328ef2a
# ╟─13b4ad4a-0599-11eb-0043-6bfc0eb23389
# ╠═e944567a-0596-11eb-2cf1-75f196bab4d1
# ╠═02673cec-0597-11eb-2db4-e9f7d9047ac0
# ╠═7e77e628-0707-11eb-1b2f-27b5c5819260
# ╠═c04a0ffe-0707-11eb-2a74-b3b258b4aec6
# ╟─358f21ca-0599-11eb-3bfb-0bb12a1465bb
# ╟─2d7ae208-0599-11eb-1a91-277b764ad45e
# ╠═c98b4f9a-0592-11eb-2034-5f05e792cbd2
# ╠═be35a06c-0596-11eb-2988-61a79c7e6a87
# ╠═a08e7230-0591-11eb-3956-fb6c2d756494
# ╟─c2640938-0591-11eb-058c-5949d1e316d3
# ╟─81629936-0593-11eb-3fca-afa68e2ace3d
# ╠═a1402cbc-0593-11eb-33b4-09a0c654eefd
# ╠═cc4209ee-0593-11eb-3762-e76539e2a455
# ╟─a398740e-0593-11eb-2f81-c92861e080b8
# ╠═dce25666-0594-11eb-24f9-85eba4949d1e
# ╟─701efa2a-0596-11eb-2592-39eced2fb875
# ╠═7c5d381a-0596-11eb-1dfe-0b7074093cd6
# ╟─0b93638a-0595-11eb-3096-af3cc7202aa5
# ╟─7c934678-0599-11eb-1e6b-1f2c9003509a
# ╠═865bbb8e-0483-11eb-0622-b9a8a6ba5d4c
# ╠═9736f388-0483-11eb-069a-0d5ed87efb73
# ╠═f39f0106-0483-11eb-3617-5fb04388489e
# ╠═52bacf14-0595-11eb-33c9-dd5c40125360
# ╠═0a5783e2-0596-11eb-1e38-a5da40325159


text:
using AtomicLevels
using WignerSymbols
using HalfIntegers
using Test

@testset "Unicode super-/subscripts" begin
    @test AtomicLevels.from_subscript("₋₊₁₂₃₄₅₆₇₈₉₀") == "-+1234567890"
    @test AtomicLevels.from_superscript("⁻⁺¹²³⁴⁵⁶⁷⁸⁹⁰") == "-+1234567890"
end

include("parity.jl")
include("orbitals.jl")
include("configurations.jl")
include("excited_configurations.jl")
include("terms.jl")
include("jj_terms.jl")
include("intermediate_terms.jl")
include("couple_terms.jl")
include("csfs.jl")
include("levels.jl")
include("jj2lsj.jl")