text:
using JLD2

# Verifies that tpf works properly with a PoolModel object
tpf_main_input = load(joinpath(dirname(@__FILE__), "../reference/tpf_poolmodel.jld2"))
tuning = tpf_main_input["tuning"]
lpd_805 = tpf_main_input["lpd_805"]
lpd_904 = tpf_main_input["lpd_904"]
data = exp.([reshape(lpd_904, 1, length(lpd_904)); reshape(lpd_805, 1, length(lpd_805))])
m = DSGE.PoolModel("ss0")

# Define state space
Φpm, Ψpm, F_ϵpm, F_upm, ~ = compute_system(m)

# Load in test inputs and outputs
test_file_inputs = load(joinpath(dirname(@__FILE__), "../reference/tpf_aux_inputs.jld2"))
test_file_outputs_pm = load(joinpath(dirname(@__FILE__), "../reference/tpf_aux_outputs_poolmodel.jld2"))

φ_old = test_file_inputs["phi_old"]
norm_weights = test_file_inputs["norm_weights"]
coeff_terms = test_file_inputs["coeff_terms"]
log_e_1_terms = test_file_inputs["log_e_1_terms"]
log_e_2_terms = test_file_inputs["log_e_2_terms"]
inc_weights = test_file_inputs["inc_weights"]
HH = ones(1,1)

s_t_nontemp = [.5; .5] * ones(1,1000) # this might be desirable since all particles same then
Random.seed!(1793)
Ψ47_pm(x) = Ψpm(x,data[:,47])
weight_kernel!(coeff_terms, log_e_1_terms, log_e_2_terms, φ_old, Ψ47_pm, data[:, 47],
               s_t_nontemp, det(HH), inv(HH);
               initialize = false, parallel = false,
               poolmodel = true)
φ_new = next_φ(φ_old, coeff_terms, log_e_1_terms, log_e_2_terms, length(data[:,47]), tuning[:r_star], 2)
correction!(inc_weights, norm_weights, φ_new, coeff_terms, log_e_1_terms, log_e_2_terms, length(data[:,47]))
true_inc_wt = φ_new^(1/2) * coeff_terms[1] * exp(log_e_1_terms[1]) * exp(φ_new * log_e_2_terms[1])

@testset "Corection and Auxiliary Tests" begin
    @test coeff_terms[1] ≈ (φ_old)^(-1)
    @test log_e_1_terms[1] ≈ log(Ψpm(s_t_nontemp[:,1], data[:,47])) * -φ_old
    @test log_e_2_terms[1] ≈ log(Ψpm(s_t_nontemp[:,1], data[:,47]))
    @test φ_new ≈ 1.0
    @test inc_weights[1] ≈ true_inc_wt
end

φ_old = test_file_inputs["phi_old"]
norm_weights = test_file_inputs["norm_weights"]
coeff_terms = test_file_inputs["coeff_terms"]
log_e_1_terms = test_file_inputs["log_e_1_terms"]
log_e_2_terms = test_file_inputs["log_e_2_terms"]
inc_weights = test_file_inputs["inc_weights"]
s_t_nontemp0 = [.5; .5] * ones(1, 900)
s_t_nontemp = hcat(s_t_nontemp0, repeat([.49; .51], 1, 100))
weight_kernel!(coeff_terms, log_e_1_terms, log_e_2_terms, φ_old, Ψ47_pm, data[:, 47],
               s_t_nontemp, det(HH), inv(HH);
               initialize = false, parallel = false,
               poolmodel = true)
φ_new = next_φ(φ_old, coeff_terms, log_e_1_terms, log_e_2_terms, length(data[:,47]), tuning[:r_star], 2)
correction!(inc_weights, norm_weights, φ_new, coeff_terms, log_e_1_terms, log_e_2_terms, length(data[:,47]))
true_inc_wt = φ_new^(1/2) * (coeff_terms[1] * exp(log_e_1_terms[1])
                             * exp(φ_new * log_e_2_terms[1]) * 900/1000
                             + coeff_terms[end] * exp(log_e_1_terms[end])
                             * exp(φ_new * log_e_2_terms[end]) * 100/1000)

@testset "Ensure different states lead to different measurement errors" begin
    @test norm_weights[1] != norm_weights[end]
    @test log_e_1_terms[1] != log_e_1_terms[end]
    @test log_e_2_terms[1] != log_e_2_terms[end]
    @test inc_weights[1] != inc_weights[end]
    @test true_inc_wt != mean(inc_weights)
end

## Mutation Tests
s_t1_temp = [.49; .51] * ones(1,1000)
ϵ_t = reshape(test_file_inputs["eps_t"][1,:], 1, 1000)
QQ = F_ϵpm.σ * ones(1,1)
accept_rate = test_file_inputs["accept_rate"]
c = test_file_inputs["c"]

c = update_c(c, accept_rate, tuning[:target_accept_rate])
Random.seed!(47)
StateSpaceRoutines.mutation!(Φpm, Ψ47_pm, QQ, det(HH), inv(HH), φ_new, data[:,47],
                             s_t_nontemp, s_t1_temp, ϵ_t, c, tuning[:n_mh_steps];
                             poolmodel = true)
@testset "Mutation Tests" begin
    @test s_t_nontemp[1] ≈ test_file_outputs_pm["s_t_nontemp"][1]
    @test ϵ_t[1] ≈ test_file_outputs_pm["eps_t"][1]
end

## Whole TPF Tests
Random.seed!(47)
s_init = reshape(rand(Uniform(.4, .6), 1000), 1, 1000)
s_init = [s_init; 1 .- s_init]
out_no_parallel = tempered_particle_filter(data, Φpm, Ψpm, F_ϵpm, F_upm, s_init;
                                           tuning..., verbose = :none, fixed_sched = [1.],
                                           parallel = false, poolmodel = true)
Random.seed!(47)
out_parallel_one_worker = tempered_particle_filter(data, Φpm, Ψpm, F_ϵpm, F_upm, s_init;
                                                   tuning..., verbose = :none, fixed_sched = [1.],
                                                   parallel = true, poolmodel = true)
@testset "TPF tests" begin
    @test out_no_parallel[1] ≈ test_file_outputs_pm["out_no_parallel"][1]
    # See tempered_particle_filter.jl's test with parallel workers
    if VERSION >= v"1.5"
        @test out_parallel_one_worker[1] ≈ -507.35270760697324
    elseif VERSION >= v"1.0"
        @test out_parallel_one_worker[1] ≈ test_file_outputs_pm["out_parallel_one_worker"][1] # should be -507.44364755284465
    end
end

nothing


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

using Markdown
using InteractiveUtils

# ╔═╡ 5bf16f26-7191-11eb-1166-73a483cf120a
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)

	using MLJ
	using CSV
	using DataFrames
	using PlutoUI
	using Test
	using Printf
	using Random
 	using Plots  # using PyPlot
end

# ╔═╡ 638e7b78-7544-11eb-2a45-9544fcbeb7bc
## include utils functions
include("./utils.jl")

# ╔═╡ 87aaecb0-7190-11eb-1572-d9110d9ffaaa
md"""
# Week 4: Ridge Regression (Interpretation)
"""

# ╔═╡ c5c92502-7190-11eb-3a75-87b7f156a7ce
md"""
In this notebook, we will run ridge regression multiple times with different L2 penalties to see which one produces the best fit. We will revisit the example of polynomial regression as a means to see the effect of L2 regularization. In particular, we will:
* Use a pre-built implementation of regression (MLJLinearModels) to run polynomial regression, this time with L2 penalty
* Use plot to visualize polynomial regressions under L2 regularization
* Choose best L2 penalty using cross-validation.
* Assess the final fit using test data.

We will continue to use the House data from previous notebooks.  (In the next programming assignment for this module, you will implement your own ridge regression learning algorithm using gradient descent.)

"""

# ╔═╡ 14d8375a-734e-11eb-2518-2d4eee6b9e3e
md"""
### Polynomial regression, revisited

We build on the material from Week 3, where we wrote the function to produce an SFrame with columns containing the powers of a given input. Copy and paste the function `polynomial_sframe` from Week 3:
"""

# ╔═╡ 55fb5a1e-734e-11eb-2019-21b125a173d1
function polynomial_df(feature; degree=3)
	@assert degree ≥ 1 "Expect degree to be ≥ 1"

	hsh = Dict{Symbol, Vector{Float64}}(:power_1 => feature)
	for deg ∈ 2:degree
		hsh[Symbol("power_$(deg)")] = feature .^ deg
	end

	return DataFrame(hsh)
end

# ╔═╡ 7a50bd00-7191-11eb-2ed6-071e6db68422
md"""
### Load in house sales data
"""

# ╔═╡ 890cb362-7191-11eb-251b-e155ed2d21f2
sales = CSV.File("../../ML_UW_Spec/C02/data/kc_house_test_data.csv";
	header=true) |> DataFrame;

# ╔═╡ b43046d2-7191-11eb-0978-b1bfe1762756
md"""
### Split data into train/test sets
"""

# ╔═╡ 76617496-734e-11eb-0177-ef53e0a61ebc
sort!(sales, [:sqft_living, :price], rev=[false, false]);

# ╔═╡ a626fab6-734e-11eb-0994-f153e8264aaf
md"""
Let us revisit the 15th-order polynomial model using the `sqft_living` input. Generate polynomial features up to degree 15 using `polynomial_df()` and fit a model with these features.

When fitting the model, use an L2 penalty of 1e-5:
"""

# ╔═╡ 9f4449e2-734e-11eb-0d61-19d14b53358a
begin
	poly_df_0 = polynomial_df(sales.sqft_living; degree=15)
	poly_df_0[!, :price] = sales.price

	first(poly_df_0, 3)
end

# ╔═╡ e985d03e-734e-11eb-26a7-effe0c92f0fc
@load RidgeRegressor pkg=MLJLinearModels

# ╔═╡ ed0179fc-734e-11eb-3309-03b0ef4cdc96
begin
	l2_small_penalty = 1e-5

	mdl0 = MLJLinearModels.RidgeRegressor(lambda=l2_small_penalty,
		penalize_intercept=false)

	X_0 = select(poly_df_0, :power_1)
	y_0 = poly_df_0.price

	mach0 = machine(mdl0, X_0, y_0)
	fit!(mach0)
	fp0 = fitted_params(mach0)

	with_terminal() do
		for (name, c) in fp0.coefs
   			println("$(rpad(name, 10)):  $(round(c, sigdigits=3))")
		end

		println("Intercept: $(round(fp0.intercept, sigdigits=3))")
	end
end

# ╔═╡ f4ac831a-7351-11eb-10b2-67efef23c436
md"""

**Quiz Qestion:  What's the learned value for the coefficient of feature `power_1`?**
  - Answer below:
"""

# ╔═╡ 4686d706-74a8-11eb-101d-c976b98872f3
## power_1
round(fp0.coefs[1][2], sigdigits=2)

# ╔═╡ 3c18cc2e-74a9-11eb-336b-1bdf4159d5f8
getfield(fp0, :coefs)[1][2]

# ╔═╡ 1fe9cdec-7352-11eb-1336-d98c99d2a002
md"""
### Observe overfitting

Recall from Week 3 that the polynomial fit of degree 15 changed wildly whenever the data changed. In particular, when we split the sales data into four subsets and fit the model of degree 15, the result came out to be very different for each subset. The model had a *high variance*. We will see in a moment that ridge regression reduces such variance.

But first, we must reproduce the experiment we did in Week 3.
"""

# ╔═╡ 4835efa8-7352-11eb-27dd-23dcf67ead55
begin
	(sales_a, sales_b) = train_test_split(sales; split=0.5, seed=42)

	(set1, set2) = train_test_split(sales_a; split=0.5, seed=42)
	(set3, set4) = train_test_split(sales_b; split=0.5, seed=42)

	(size(set1), size(set2), size(set3), size(set4))
end

# ╔═╡ b3b02ed6-7352-11eb-3a15-9f392deb56d8
begin

function make_poly(tset; degree=15, output=:price)
  poly_df = polynomial_df(tset.sqft_living; degree)
  features = names(poly_df)
  poly_df[!, output] = tset[!, output]
  (features, poly_df)
end

function fit_poly(tset; degree=15, output=:price, l2_penalty=l2_small_penalty)
  (features, poly_df) = make_poly(tset; degree, output)
  mdl = MLJLinearModels.RidgeRegressor(lambda=l2_penalty,
		penalize_intercept=false)
  X_ = select(poly_df, features)
  y_ = poly_df[!, output]

  mach = machine(mdl, X_, y_)
  fit!(mach)

  (mach, X_, y_)
end

function print_coeff(mach, ix)
  fp = fitted_params(mach)
  println("\n- model $(ix)")
  for (name, c) in fp.coefs[1:1]  # print only first coeff
	println("$(rpad(name, 10)):  $(round(c, sigdigits=3))")
  end

  println("Intercept: $(round(fp.intercept, sigdigits=3))")
end

function find(machs; key=:min)
	fps = map(m -> fitted_params(m), machs) |>
		fps -> map(fp -> getfield(fp, :coefs)[1][2], fps)

	ix = key == :min ? argmin(fps) : argmax(fps)
	(fps[ix], string("model_", ix))
end

end

# ╔═╡ 6431250e-7354-11eb-1446-37f9f51cb959
md"""
##### Set1, Set2, Set3 & Set4
"""

# ╔═╡ 64f0b460-7353-11eb-23b5-0b2a0c9006f4
begin
	hsh = Dict{Symbol, Vector{Any}}(:machine => [], :tset => [])

	for (ix, s) ∈ enumerate((set1, set2, set3, set4))
		(mach_, X_, y_) = fit_poly(s)
		push!(hsh[:machine], mach_)
		push!(hsh[:tset], (X_, y_))
	end

	ps = []

	for ix ∈ 1:length(hsh[:machine])
		(X_, y_) = hsh[:tset][ix]
		mach_ = hsh[:machine][ix]

		p1 = scatter(X_.power_1, y_, legend=false, color=[:lightblue], marker=".")
 		p2 = scatter!(X_.power_1, predict(mach_, X_), color=[:orange], marker="-")
		push!(ps, (p1, p2))
	end

end

# ╔═╡ c7c9f6dc-7495-11eb-2bd5-8ba7f20f997a
with_terminal() do
	for ix in 1:length(hsh[:machine])
		print_coeff(hsh[:machine][ix], ix)
	end
end

# ╔═╡ e22363b4-7499-11eb-3b50-5905e5eef045
begin
	lg1 = grid(2, 2, widths=[0.0, 0.9], heights=[0.5, 0.5])
	plot(ps[1]..., ps[2]..., layout=lg1)
end

# ╔═╡ f08611a8-749a-11eb-1014-e9e29e34a966
begin
	lg2 = grid(2, 2, widths=[0.0, 0.9], heights=[0.5, 0.5])
	plot(ps[3]..., ps[4]..., layout=lg2)
end

# ╔═╡ e4c6e384-7354-11eb-307d-579569191332
md"""
The four curves should differ from one another a lot, as should the coefficients you learned.

**Quiz Question:  For the models learned in each of these training sets, what are the smallest and largest values you learned for the coefficient of feature `power_1`?**

  - cf. below
"""

# ╔═╡ ac0d2d96-74a8-11eb-079e-036266a3b75e
(smallest=find(hsh[:machine]), largest=find(hsh[:machine]; key=:max))

# ╔═╡ 40bd2ee6-7355-11eb-00ba-1943adaa3fb9
md"""
### Ridge regression comes to rescue

Generally, whenever we see weights change so much in response to change in data, we believe the variance of our estimate to be large. Ridge regression aims to address this issue by penalizing "large" weights. (Weights of `model15` looked quite small, but they are not that small because 'sqft_living' input is in the order of thousands.)

With the argument `l2_penalty=1e5`, fit a 15th-order polynomial model on `set1`, `set2`, `set3`, and `set4`. Other than the change in the `l2_penalty` parameter, the code should be the same as the experiment above.
"""

# ╔═╡ 5aa44ab0-7355-11eb-17cb-7bd692a55a37
md"""
#### Set1, Set2, Set3 & Set4
"""

# ╔═╡ 5aa4d1cc-7355-11eb-09a9-dd185f7871e1
l2_penalty=1e5

# ╔═╡ 96c23cd0-749c-11eb-0c80-63c4f283c2e5
begin
	hsh_v1r = Dict{Symbol, Vector{Any}}(:machine => [], :tset => [])

	for (ix, s) ∈ enumerate((set1, set2, set3, set4))
		(mach_, X_, y_) = fit_poly(s; l2_penalty)
		push!(hsh_v1r[:machine], mach_)
		push!(hsh_v1r[:tset], (X_, y_))
	end

	## plot prep.
	ps2 = []
	for ix ∈ 1:length(hsh[:machine])
		(X_, y_) = hsh[:tset][ix]
		mach_ = hsh[:machine][ix]

		p1 = scatter(X_.power_1, y_, legend=false, color=[:lightblue], marker=".")
 		p2 = scatter!(X_.power_1, predict(mach_, X_), color=[:orange], marker="-")
		push!(ps2, (p1, p2)) # plot(p1, p2, layout=l)
	end

	with_terminal() do
		for ix in 1:length(hsh_v1r[:machine])
			print_coeff(hsh_v1r[:machine][ix], ix)
		end
	end
end

# ╔═╡ 03e1ed82-74a8-11eb-140c-e31183a2314b
begin
	lg11 = grid(2, 2, widths=[0.0, 0.9], heights=[0.5, 0.5])
	plot(ps2[1]..., ps2[2]..., layout=lg11)
end

# ╔═╡ 034e5e96-74a8-11eb-030d-fb9f7ea9af55
begin
	lg12 = grid(2, 2, widths=[0.0, 0.9], heights=[0.5, 0.5])
	plot(ps2[3]..., ps2[4]..., layout=lg12)
end

# ╔═╡ d32b7d62-7355-11eb-00f5-e15ead474ce5
md"""
These curves should vary a lot less, now that you applied a high degree of regularization.

**Quiz Question:  For the models learned with the high level of regularization in each of these training sets, what are the smallest and largest values you learned for the coefficient of feature `power_1`?**

  - cf. below
"""

# ╔═╡ a96763a2-74aa-11eb-2a92-87ee73d05221
(smallest=find(hsh_v1r[:machine]), largest=find(hsh_v1r[:machine]; key=:max))

# ╔═╡ 8e9ed252-7357-11eb-117c-91f7ff585a2d
md"""
### Selecting an L2 penalty via cross-validation

Just like the polynomial degree, the L2 penalty is a "magic" parameter we need to select. We could use the validation set approach as we did in the last module, but that approach has a major disadvantage: it leaves fewer observations available for training. **Cross-validation** seeks to overcome this issue by using all of the training set in a smart way.

We will implement a kind of cross-validation called **k-fold cross-validation**. The method gets its name because it involves dividing the training set into k segments of roughtly equal size. Similar to the validation set method, we measure the validation error with one of the segments designated as the validation set. The major difference is that we repeat the process k times as follows:

Set aside segment 0 as the validation set, and fit a model on rest of data, and evalutate it on this validation set

Set aside segment 1 as the validation set, and fit a model on rest of data, and evalutate it on this validation set

...

Set aside segment k-1 as the validation set, and fit a model on rest of data, and evalutate it on this validation set

After this process, we compute the average of the k validation errors, and use it as an estimate of the generalization error. Notice that  all observations are used for both training and validation, as we iterate over segments of data.

To estimate the generalization error well, it is crucial to shuffle the training data before dividing them into segments.

"""

# ╔═╡ 91a058b6-7359-11eb-1da5-eba5d339d100
function get_rss(mach, X, y)
    ŷ = predict(mach, X)     # First get the predictions
    diff = y .- ŷ            # Then compute the residuals/errors
    rss = sum(diff .* diff)  # Then square and add them up
    return rss
end

# ╔═╡ 957b90fe-7359-11eb-23fc-039055cd2914
function k_fold_cross_validation(k, l2_penalty, df, output, features;
	degree=15)
	n = size(df)[1]
	rss = []
	for ix ∈ 1:k
		start_ = (n * (ix - 1)) ÷ k + 1
		end_ = (n * ix) ÷ k

		## get the partition train/valid
		valid_set = df[start_:end_, :]
		train_set = vcat(df[1:start_ - 1, :], df[end_+1:end, :])

		##
		mdl = MLJLinearModels.RidgeRegressor(lambda=l2_penalty,
			penalize_intercept=false)
  		X_ = select(train_set, features)
  		y_ = train_set[!, output]

  		mach = machine(mdl, X_, y_)
  		fit!(mach)

		##
		push!(rss, get_rss(mach, valid_set, valid_set[!, output]))
	end
	sum(rss) / k
end

# ╔═╡ 9af5697a-73af-11eb-1558-dbc635b07638
function find_best_l2(train_valid; k=10, degree=15, output=:price)
	l2_rss, l2_p = [], []
	best_l2, min_rss = nothing, nothing
	(features, poly_set) = make_poly(train_valid; degree, output)

	for cur_l2p ∈ (10^ix for ix ∈ 1:0.5:7)
		push!(l2_p, cur_l2p)
		rss = k_fold_cross_validation(k, cur_l2p, poly_set, output, features)
		push!(l2_rss, rss)

		# @printf("current rss: %2.5e / current penalty: %2.5f\n", rss, cur_l2p)
		if isnothing(min_rss) || rss < min_rss
			min_rss = rss
			best_l2 = cur_l2p
		end
	end
	(best_l2, min_rss, l2_p, l2_rss)
end

# ╔═╡ 51150d2e-7358-11eb-2d17-c5c6898b7045
begin
	# train_valid is shuffled during the train-test_split
	(train_valid, test) = train_test_split(sales; split=0.9, seed=42)

	k = 10
	degree = 15
	output=:price

	(best_l2, min_rss, l2_p, l2_rss) = find_best_l2(train_valid; k, degree, output)
	with_terminal() do
		@printf("min rss: %2.5e / best penalty: %2.5e\n", min_rss, best_l2)
	end
end

# ╔═╡ daa7efec-73b0-11eb-3657-d520436222d5
plot(l2_p, l2_rss, marker="-", color=[:lightblue], xscale=:log)

# ╔═╡ b273d2d6-73b1-11eb-06df-1559778e561d
md"""
Once you found the best value for the L2 penalty using cross-validation, it is important to retrain a final model on all of the training data using this value of l2_penalty.

This way, your final model will be trained on the entire dataset.
"""

# ╔═╡ 199f8240-73b2-11eb-2753-ab74e806ab9f
begin
	(features, train_set) = make_poly(train_valid; degree, output)

	mdl = MLJLinearModels.RidgeRegressor(lambda=best_l2,
			penalize_intercept=false)
  		X_ = select(train_set, features)
  		y_ = train_set[!, output]

  		mach = machine(mdl, X_, y_)
  		fit!(mach)

	(_, test_set) = make_poly(test; degree, output)

	final_rss = get_rss(mach, test_set, test_set[!, output])

	with_terminal() do
		@printf("final rss: %2.5e\n", final_rss)
	end
end

# ╔═╡ 3db30566-73b3-11eb-1d4e-f919ad63641a
md"""
**Quiz Question: Using the best L2 penalty found above, train a model using all *training* data. What is the RSS on the *test* data of the model you learn with this L2 penalty?**

 - Final rss: 4.03826e+14
"""

# ╔═╡ Cell order:
# ╟─87aaecb0-7190-11eb-1572-d9110d9ffaaa
# ╟─c5c92502-7190-11eb-3a75-87b7f156a7ce
# ╠═5bf16f26-7191-11eb-1166-73a483cf120a
# ╟─14d8375a-734e-11eb-2518-2d4eee6b9e3e
# ╠═55fb5a1e-734e-11eb-2019-21b125a173d1
# ╟─7a50bd00-7191-11eb-2ed6-071e6db68422
# ╠═890cb362-7191-11eb-251b-e155ed2d21f2
# ╠═638e7b78-7544-11eb-2a45-9544fcbeb7bc
# ╟─b43046d2-7191-11eb-0978-b1bfe1762756
# ╠═76617496-734e-11eb-0177-ef53e0a61ebc
# ╟─a626fab6-734e-11eb-0994-f153e8264aaf
# ╠═9f4449e2-734e-11eb-0d61-19d14b53358a
# ╠═e985d03e-734e-11eb-26a7-effe0c92f0fc
# ╠═ed0179fc-734e-11eb-3309-03b0ef4cdc96
# ╟─f4ac831a-7351-11eb-10b2-67efef23c436
# ╠═4686d706-74a8-11eb-101d-c976b98872f3
# ╠═3c18cc2e-74a9-11eb-336b-1bdf4159d5f8
# ╟─1fe9cdec-7352-11eb-1336-d98c99d2a002
# ╠═4835efa8-7352-11eb-27dd-23dcf67ead55
# ╠═b3b02ed6-7352-11eb-3a15-9f392deb56d8
# ╟─6431250e-7354-11eb-1446-37f9f51cb959
# ╠═64f0b460-7353-11eb-23b5-0b2a0c9006f4
# ╠═c7c9f6dc-7495-11eb-2bd5-8ba7f20f997a
# ╠═e22363b4-7499-11eb-3b50-5905e5eef045
# ╠═f08611a8-749a-11eb-1014-e9e29e34a966
# ╟─e4c6e384-7354-11eb-307d-579569191332
# ╠═ac0d2d96-74a8-11eb-079e-036266a3b75e
# ╟─40bd2ee6-7355-11eb-00ba-1943adaa3fb9
# ╟─5aa44ab0-7355-11eb-17cb-7bd692a55a37
# ╠═5aa4d1cc-7355-11eb-09a9-dd185f7871e1
# ╠═96c23cd0-749c-11eb-0c80-63c4f283c2e5
# ╠═03e1ed82-74a8-11eb-140c-e31183a2314b
# ╠═034e5e96-74a8-11eb-030d-fb9f7ea9af55
# ╟─d32b7d62-7355-11eb-00f5-e15ead474ce5
# ╠═a96763a2-74aa-11eb-2a92-87ee73d05221
# ╟─8e9ed252-7357-11eb-117c-91f7ff585a2d
# ╠═91a058b6-7359-11eb-1da5-eba5d339d100
# ╠═957b90fe-7359-11eb-23fc-039055cd2914
# ╠═9af5697a-73af-11eb-1558-dbc635b07638
# ╠═51150d2e-7358-11eb-2d17-c5c6898b7045
# ╠═daa7efec-73b0-11eb-3657-d520436222d5
# ╟─b273d2d6-73b1-11eb-06df-1559778e561d
# ╠═199f8240-73b2-11eb-2753-ab74e806ab9f
# ╟─3db30566-73b3-11eb-1d4e-f919ad63641a

ext: jl


text:
#=
My very basic visualization toolkit for the VEM using the Plots module.
=#

using Plots
pyplot() # select PyPlot background

#= The next three functions are used to produce the Jet colormap. =#
function interpolate(v::Float64, y0::Float64, x0::Float64, y1::Float64, x1::Float64)
    return (v - x0)*(y1 - y0)/(x1 - x0) + y0
end

function base(v::Float64)
    if v <= -0.75
        return 0.
    elseif v <= -0.25
        return interpolate(v, 0., -0.75, 1., -0.25)
    elseif v <= 0.25
        return 1.
    elseif v <= 0.75
        return interpolate(v, 1., 0.25, 0., 0.75)
    else
        return 0.
    end
end

red(gs::Float64)   = base(gs - 0.5)
green(gs::Float64) = base(gs)
blue(gs::Float64)  = base(gs + 0.5)

function plotSolutionNodal(msh::mesh, u::Array{Float64, 1}; myTitle::String="",
    toFile::Bool=false)
    #= Plots nodal values =#
    println("Printing solution")
    # set polygons
    nvtx = size(msh.vtx, 2)
    nelt = length(msh.elt)
    moyU = zeros(Float64, nelt)
    # set polygons
    el = msh.elt[1]
    xl = collect(transpose(msh.vtx[:,el]))
    println
    pols = [xl; NaN NaN]
    ul   = u[el]
    moyU[1] = sum(ul)/length(ul)
    for i=2:nelt
        @inbounds el      = msh.elt[i]
        @inbounds xl      = collect(transpose(msh.vtx[:, el]))
        tmp = [xl; NaN NaN]
        pols = [pols; tmp]
        @inbounds ul      = u[el]
        @inbounds moyU[i] = sum(ul)/length(ul)
    end
    # set colors
    maxU = maximum(moyU)
    minU = minimum(moyU)
    delU = maxU - minU
    gs   = -1 + 2/delU*(moyU[1] - minU)
    colU = [RGBA(red(gs), green(gs), blue(gs), 1)]
    for i=2:nelt
        @inbounds gs = -1 + 2/delU*(moyU[i] - minU)
        push!(colU, RGBA(red(gs), green(gs), blue(gs), 1))
    end
    # plot
    fig = Plots.plot(pols[:,1], pols[:,2], seriestype=:shape, c=colU, show=true,
    label=nothing, linewidth=0.2, linecolor=RGBA(255/255,255/255,255/255,1),
    aspect_ratio=:equal, title=myTitle)
    if toFile
        png(fig, "myVEM")
    end
    display(fig)
    nothing
end

function plotSolutionElt(msh::mesh, ul::Array{Float64, 1}; myTitle::String="",
    toFile::Bool=false)
    #= Plots a function constant per element, for example the pressure field. =#
    println("Printing solution")
    # set polygons
    nvtx = size(msh.vtx, 2)
    nelt = length(msh.elt)
    moyU = zeros(Float64, nelt)
    # set polygons
    el = msh.elt[1]
    xl = collect(transpose(msh.vtx[:,el]))
    println
    pols = [xl; NaN NaN]
    for i=2:nelt
        @inbounds el      = msh.elt[i]
        @inbounds xl      = collect(transpose(msh.vtx[:, el]))
        tmp = [xl; NaN NaN]
        pols = [pols; tmp]
    end
    # set colors
    maxU = maximum(ul)
    minU = minimum(ul)
    delU = maxU - minU
    gs   = -1 + 2/delU*(ul[1] - minU)
    colU = [RGBA(red(gs), green(gs), blue(gs), 1)]
    for i=2:nelt
        @inbounds gs = -1 + 2/delU*(ul[i] - minU)
        push!(colU, RGBA(red(gs), green(gs), blue(gs), 1))
    end
    #
    fig = Plots.plot(pols[:,1], pols[:,2], seriestype=:shape, c=colU, show=true,
    label=nothing, linewidth=0.2, linecolor=RGBA(255/255,255/255,255/255,1),
    aspect_ratio=:equal, title=myTitle)
    if toFile
        png(fig, "myVEM")
    end
    display(fig)
    nothing
end


text:
module TestExamples1DMHD

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1d_dgsem")

@testset "MHD Multicomponent" begin

  @trixi_testset "elixir_mhdmulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [0.08160481582829862, 0.05467911944326103, 0.05467911944326103, 0.15450926504459692,
              4.130462730494e-17, 0.054763752050210085, 0.054763752050210085, 0.008371564857135208,
              0.016743129714270416, 0.03348625942854083],
      linf = [0.18198258075330706, 0.09136114386311774, 0.09136114386311774, 0.423831369951313,
              1.1102230246251565e-16, 0.09937317613143604, 0.09937317613143604, 0.0157164284712992,
              0.0314328569425984, 0.0628657138851968])
  end

  @trixi_testset "elixir_mhdmulti_ec.jl with flux_derigs_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [0.08153372259925547, 0.05464109003345891, 0.05464109003345891, 0.1540576724164453,
              4.130462730494e-17, 0.054734930802131036, 0.054734930802131036, 0.008391254781284321,
              0.016782509562568642, 0.033565019125137284],
      linf = [0.17492544007323832, 0.09029632168248182, 0.09029632168248182, 0.40798609353896564,
              1.1102230246251565e-16, 0.09872923637833075, 0.09872923637833075, 0.01609818847160674,
              0.03219637694321348, 0.06439275388642696],
      volume_flux = flux_derigs_etal)
  end

  @trixi_testset "elixir_mhdmulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"),
      l2   = [0.07968782477167513, 0.05398115008116676, 0.05398115008116676, 0.15015281822439228,
              4.130462730494e-17, 0.053629890024921495, 0.053629890024921495, 0.008279068245579706,
              0.016558136491159413, 0.033116272982318826],
      linf = [0.14118014632124837, 0.07820697032983395, 0.07820697032983395, 0.3390558674728652,
              1.1102230246251565e-16, 0.06998787893467828, 0.06998787893467828, 0.014943825414763745,
              0.02988765082952749, 0.05977530165905498])
  end

  @trixi_testset "elixir_mhdmulti_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_convergence.jl"),
      l2   = [1.7337265267786785e-5, 0.00032976971029271364, 0.0003297697102926479, 6.194071694759044e-5,
              4.130462730494001e-17, 0.00032596825025664136, 0.0003259682502567132, 2.5467510126885455e-5,
              5.093502025377091e-5, 0.00010187004050754182],
      linf = [3.877554303711845e-5, 0.0012437848638874956, 0.0012437848638876898, 0.00016431262020277781,
              1.1102230246251565e-16, 0.0012443734922607112, 0.001244373492260704, 5.691007974162332e-5,
              0.00011382015948324664, 0.00022764031896649328])
    end

  @trixi_testset "elixir_mhdmulti_briowu_shock_tube.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_briowu_shock_tube.jl"),
      l2   = [0.1946577804333822, 0.3591196215528672, 0.0, 0.36875476066849383,
              4.7644020131827105e-16, 0.36668249926193885, 0.0, 0.05775369214541893,
              0.11550738429083786],
      linf = [0.4345551123140612, 1.0874941615375844, 0.0, 1.0493729052116585,
              3.219646771412954e-15, 1.5160434573973656, 0.0, 0.18616213071936066,
              0.3723242614387213])
    end

end

end # module


text:
# 2-arg show: show schema and list # of metadata entries if non-zero
function Base.show(io::IO, table::Table)
    ncols = length(Tables.columnnames(table))
    print(io, "$(typeof(table)) with $(Tables.rowcount(table)) rows, $(ncols) columns,")
    meta = getmetadata(table)
    if meta !== nothing && !isempty(meta)
        print(io, " ", length(meta), " metadata entries,")
    end
    sch = Tables.schema(table)
    print(io, " and schema:\n")
    show(IOContext(io, :print_schema_header => false), sch)
    return nothing
end

# 3-arg show: show schema and show metadata entries adaptively according to `displaysize`
function Base.show(io::IO, mime::MIME"text/plain", table::Table)
    limit = get(io, :limit, false)::Bool
    display_rows, display_cols = displaysize(io)
    ncols = length(Tables.columnnames(table))
    meta = getmetadata(table)
    if meta !== nothing
        display_rows -= 1 # decrement for metadata header line
        display_rows -= min(length(meta), 2) # decrement so we can show at least 2 lines of metadata
    end
    print(io, "$(typeof(table)) with $(Tables.rowcount(table)) rows, $(ncols) columns, and ")
    sch = Tables.schema(table)
    print(io, "schema:\n")
    schema_context = IOContext(io, :print_schema_header => false, :displaysize => (max(display_rows, 3), display_cols))
    schema_str = sprint(show, mime, sch; context=schema_context)
    print(io, schema_str)
    display_rows -= (count("\n", schema_str) + 1) # decrement for number of lines printed
    if meta !== nothing
        print(io, "\n\nwith metadata given by a ")
        show(IOContext(io, :displaysize => (max(display_rows, 5), display_cols)), mime, meta)
    end
    return nothing
end


text:
module Monkey

include("lexer/lexer.jl")
include("parser/parser.jl")

export Lexer, Parser, parse

end

max_stars_repo_name: viert/monkeylang-julia

max_issues_repo_name: viert/monkeylang-julia

max_forks_repo_name: viert/monkeylang-julia


text:
module AnalyticHierarchyProcess

using Compat:view

export eigenv,CRA,lmax,ctriads,entropia,vikor

# Funtion lambda max###############################################################
function lmax(x)
data = x
Z=data
A=sum(data,1)
Anorm=data./A
tam=length(Z[1,:])
wprom=sum(Anorm,2)/tam
Awprom=Z*wprom
lambdamax=mean(Awprom./wprom)
return lambdamax
end

#Function weigth###############################################################
function eigenv(x)
data = x
Z=data
A=sum(data,1)
Anorm=data./A
tam=length(Z[1,:])
wprom=sum(Anorm,2)/tam
Awprom=Z*wprom
lambdamax=mean(Awprom./wprom)

if tam==1
                      CI=0
           elseif tam==2
                      CI=0
           else
           CI=(lambdamax-tam)/(tam-1)
           end

if tam==1
    CR=0.00
elseif tam==2
    CR=0.00
elseif tam==3
    CR=CI/0.52
elseif tam==4
    CR=CI/0.89
elseif tam==5
    CR=CI/1.11
elseif tam==6
    CR=CI/1.25
elseif tam==7
    CR=CI/1.35
elseif tam==8
    CR=CI/1.40
elseif tam==9
    CR=CI/1.45
elseif tam==10
    CR=CI/1.49
end

if tam>=3
M=Z^100
sumafilas=sum(M,2)
WCRIT=sumafilas/sum(sumafilas)
else
M=Z
sumafilas=sum(M,2)
WCRIT=sumafilas/sum(sumafilas)

end

return   WCRIT

end

function CRA(x)  ############################################################################################
data = x
Z=data
A=sum(data,1)
Anorm=data./A
tam=length(Z[1,:])
wprom=sum(Anorm,2)/tam
Awprom=Z*wprom
lambdamax=mean(Awprom./wprom)

if tam==1
                      CI=0
           elseif tam==2
                      CI=0
           else
           CI=(lambdamax-tam)/(tam-1)
           end

if tam==1
    CR=0.00
elseif tam==2
    CR=0.00
elseif tam==3
    CR=CI/0.52
elseif tam==4
    CR=CI/0.89
elseif tam==5
    CR=CI/1.11
elseif tam==6
    CR=CI/1.25
elseif tam==7
    CR=CI/1.35
elseif tam==8
    CR=CI/1.40
elseif tam==9
    CR=CI/1.45
elseif tam==10
    CR=CI/1.49
end

return CR
end
function ctriads(t,Criteria) #####################################################################
     c=[]
for T=1:length(t)

ai=Criteria[t[T][1], t[T][2]]
aj=Criteria[t[T][1], t[T][3]]
ak=Criteria[t[T][2], t[T][3]]

m=[1    ai   aj
  1/ai  1    ak
  1/aj  1/ak 1]
 co=CRA(m)
push!(c,co)
end

return c

end
################################################
function entropia(m)
a=size(m,1)
    k=-1/log(a)
    mnorm=m./sum(m,1)
    mlog=mnorm.*log.(mnorm)
    sumcol=sum(mlog,1)
    entropia1=sumcol*k
    diversidad=1-entropia1
    w=(diversidad/sum(diversidad,2))
     return w
end
#################################################
function vikor(a,wi,minmax)
tamano=size(a,2) #saber cuantos criterios se tienen
    filas=size(a,1)
    vaciomejor=zeros(1,tamano)
    vaciopeor=zeros(1,tamano)

    for j=1:tamano
        col=a[:,j]
        if minmax[j]==1
            vaciomejor[j]=maximum(col)
            vaciopeor[j]=minimum(col)
        elseif minmax[j]==0
            vaciomejor[j]=minimum(col)
            vaciopeor[j]=maximum(col)
        end
    end
    println("El vector mejor es=",vaciomejor)
    println("El vector peor es=",vaciopeor)

    matrizrs=zeros(filas,tamano)
    for j=1:tamano
        for i=1:filas
            matrizrs[i,j]=
            wi[j]*((vaciomejor[j]-a[i,j])/(vaciomejor[j]-vaciopeor[j]))
        end
    end
    println("La matriz rs es=",matrizrs)
    sj=sum(matrizrs,2)
    rj=maximum(matrizrs,2)
    println("El vector sj=",sj)
    println("El vector rj=",rj)
    qj=zeros(filas,1)
    for i=1:filas
        qj[i]=(0.5*(sj[i]-minimum(sj))/(maximum(sj)-minimum(sj)))+
        (0.5*(rj[i]-minimum(rj))/(maximum(rj)-minimum(rj)))
    end
    println("el indice qj=",qj)
return qj
end
end # module


text:
#
# immutable RUsage
#     ru_utime_sec::Clong         #  user CPU time used
#     ru_utime_usec::Clong        #  user CPU time used
#     ru_stime_sec::Clong         #  system CPU time used
#     ru_stime_usec::Clong        #  system CPU time used
#     ru_maxrss::Clong            #  maximum resident set size
#     ru_ixrss::Clong             #  integral shared memory size
#     ru_idrss::Clong             #  integral unshared data size
#     ru_isrss::Clong             #  integral unshared stack size
#     ru_minflt::Clong            #  page reclaims (soft page faults)
#     ru_majflt::Clong            #  page faults (hard page faults)
#     ru_nswap::Clong             #  swaps
#     ru_inblock::Clong           #  block input operations
#     ru_oublock::Clong           #  block output operations
#     ru_msgsnd::Clong            #  IPC messages sent
#     ru_msgrcv::Clong            #  IPC messages received
#     ru_nsignals::Clong          #  signals received
#     ru_nvcsw::Clong             #  voluntary context switches
#     ru_nivcsw::Clong            #  involuntary context switches
# end
#
# RUSAGE_SELF::Cint is defined a 0 both on Mac and Linux....
#
# ru1 = Array(RUsage, 1)
# ru2 = Array(RUsage, 1)
#
# diff = ru2[1].ru_maxrss - ru1[1].ru_maxrss
# WARN = (diff > 1000) ? "WARNING" : ""
#
# ccall(:getrusage, Cint, (Cint, Ptr{Void}), 0, ru1)
# ccall(:getrusage, Cint, (Cint, Ptr{Void}), 0, ru2)

function get_vmsize()
    lines = split(open(readall, "/proc/self/status"), "\n")
    vm_line = split(filter(x->beginswith(x, "VmSize"), lines)[1])
    vmsize = int(vm_line[2])

    if uppercase(vm_line[3]) == "KB"
        vmunits = 1
    elseif uppercase(vm_line[3]) == "MB"
        vmunits = 1000
    else
        error("Unknown unit $(vm_line[3])")
    end

    vmsize * vmunits
end

function run_mtest(name, testf)
    print ("Testing $name...")
    for i in 1:2
        print ("priming process...")
        testf()             # priming
    end
    vm1 = get_vmsize()
    println ("monitored run...")
    testf()             # actual run
    vm2 = get_vmsize()

    diff = vm2 - vm1
    WARN = (diff > 1000) ? "<===================================== WARNING" : ""
    println("Memory Test ($name) : VMSize difference : $diff KB $WARN")

end

function mtest_create_strings()
    for i in 1:10^8
        string("$i")
    end

    gc()
end

run_mtest("create_strings", () -> mtest_create_strings())

function mtest_remotecall_fetch()
    for i in 1:10^5
        remotecall_fetch(1, myid)
    end
    gc()
end

run_mtest("remotecall_fetch", () -> mtest_remotecall_fetch())

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function window_signals!(widgets, signals)
    signals[:slider] = map(GtkReactive.value, widgets[:slider])

    button_signal = GtkReactive.Signal(false)
    Gtk.signal_connect(widgets[:button].widget, :clicked) do widget
        push!(button_signal, true)
    end
    signals[:button_done] = button_signal
end

function create_window(name::String, canvassize::Tuple{Int,Int})
    gridsize = (1, 1)
    winsize  = ImageView.canvas_size(ImageView.screen_size() ./ 2,
                                     map(*, canvassize, gridsize))
    win     = Gtk.Window(name, winsize...)
    ImageView.window_wrefs[win] = nothing
    Gtk.signal_connect(win, :destroy) do w
        exit()
    end

    vbox = Gtk.Box(:v)
    push!(win, vbox)

    menubar, menu_signals = create_menubar(win)
    push!(vbox, menubar)

    frames, canvases = ImageView.frame_canvas(:auto)
    push!(vbox, frames)

    video_range   = 1:100
    player_signal = GtkReactive.Signal(1)
    player_widget = GtkReactive.player(player_signal, video_range)

    grid = Gtk.Grid()
    Gtk.set_gtk_property!(grid, :column_homogeneous, true)
    slider_widget = GtkReactive.slider(0.0:0.01:1.0)
    button_widget = GtkReactive.button("done")
    grid[1:4,1] = player_widget
    grid[1:3,2] = slider_widget
    grid[4  ,2] = button_widget

    push!(vbox, grid)

    widgets = Dict(:window=>win,
                   :vbox=>vbox,
                   :grid=>grid,
                   :frame=>frames,
                   :slider=>slider_widget,
                   :button=>button_widget,
                   :player=>player_widget,
                   :canvas=>canvases,
                   :player_range=>video_range)
    window_signals = Dict{Symbol, GtkReactive.Signal}(:player=>player_signal)
    signals        = merge(window_signals, menu_signals)

    window_signals!(widgets, signals)

    widgets, signals
end


text:
module Sugars
using MLStyle
using MLStyle.AbstractPatterns

export Q, And, Or
@nospecialize

struct Q end
struct And end
struct Or end

function MLStyle.pattern_unref(
    ::Type{Or},
    self::Function,
    args::AbstractArray
)
    isempty(args) && error("An Or pattern should take more than 1 clause!")
    or([self(arg) for arg in args])
end

function MLStyle.pattern_unref(
    ::Type{And},
    self::Function,
    args::AbstractArray
)
    isempty(args) && error("An And pattern should take more than 1 clause!")
    and([self(arg) for arg in args])
end

function MLStyle.pattern_unref(
    ::Type{Q},
    self::Function,
    args::AbstractArray
)
    length(args) === 1 || error("A Q pattern should take only 1 argument!")
    self(Expr(:quote, args[1]))
end
@specialize
end


text:
"""
    cor_nearPD(R::Matrix{T}, τ::Real=1e-6; tol::Real=1e-3) where {T<:AbstractFloat}

Return the nearest positive definite correlation matrix to `R`. `τ` is a
tuning parameter that controls the minimum eigenvalue of the resulting matrix.
`τ` can be set to zero if only a positive semidefinite matrix is needed. `tol`
is the accuracy in the relative gap. Set to a smaller value (e.g. `1e-8`) if
higher accuracy is needed.

See also: [`cor_fastPD`](@ref), [`cor_fastPD!`](@ref)

# Examples
```jldoctest
julia> import LinearAlgebra: isposdef

julia> r = [1.00 0.82 0.56 0.44; 0.82 1.00 0.28 0.85; 0.56 0.28 1.00 0.22; 0.44 0.85 0.22 1.00]
4×4 Array{Float64,2}:
 1.0   0.82  0.56  0.44
 0.82  1.0   0.28  0.85
 0.56  0.28  1.0   0.22
 0.44  0.85  0.22  1.0

julia> isposdef(r)
false

julia> p = cor_nearPD(r)
4×4 Array{Float64,2}:
 1.0       0.817095  0.559306  0.440514
 0.817095  1.0       0.280196  0.847352
 0.559306  0.280196  1.0       0.219582
 0.440514  0.847352  0.219582  1.0

julia> isposdef(p)
true
```
"""
function cor_nearPD(R::Matrix{T}, τ::Real=1e-6; tol::Real=1e-3) where {T<:AbstractFloat}

    # Setup
    n = size(R, 1)
    iter_outer = 200
    iter_inner = 20
    iter_cg    = 200
    tol_cg     = 1e-2
    tol_ls     = 1e-4
    err_tol    = max(eps(T), tol)

    # Make R symmetric
    R .= Symmetric(R, :U) # [n,n]
    R[diagind(R)] .= one(T)

    b = ones(T, n)
    if τ > 0
        b .-= τ
        R[diagind(R)] .-= τ
    end
    b₀ = copy(b)

    y    = zeros(T, n)  # [n,1]
    X    = copy(R)            # [n,n]
    λ, P = eigen(X)           # [n,1], [n,n]
    λ, P = Vector{T}(λ), Matrix{T}(P)
    λ   .= reverse(λ)         # [n,1]
    P   .= reverse(P, dims=2) # [n,n]

    f₀, Fy = npd_gradient(y, λ, P, b₀, n) # [1], [n,1]
    f      = f₀      # [1]
    b     .= b₀ - Fy # [n,1]

    Ω₀ = npd_set_omega(λ, n) # [r,s]
    x₀ = copy(y)             # [n,1]

    X       .= npd_pca(b₀, X, λ, P, n) # [n,n]
    val_R    = 0.5 * norm(R)^2
    val_dual = val_R - f₀
    val_obj  = 0.5 * norm(X - R)^2
    gap      = (val_obj - val_dual) / (1 + abs(val_dual) + abs(val_obj))

    norm_b  = norm(b)
    norm_b0 = norm(b₀) + 1
    norm_b_rel = norm_b / norm_b0

    k = 0
    while (gap > err_tol) && (norm_b_rel > err_tol) && (k < iter_outer)
        c = npd_precond_matrix(Ω₀, P, n)                # [n,1]
        d = npd_pre_cg(b, c, Ω₀, P, tol_cg, iter_cg, n) # [n,1]

        slope = sum((Fy - b₀) .* d)          # [1]
        y    .= x₀ + d                       # [n,1]
        X    .= R + diagm(y)                 # [n,n]
        λ, P  = eigen(X)                     # [n,1], [n,n]
        λ, P = Vector{T}(λ), Matrix{T}(P)
        λ    .= reverse(λ)                   # [n,1]
        P    .= reverse(P, dims=2)           # [n,n]
        f, Fy = npd_gradient(y, λ, P, b₀, n) # [1], [n,1]

        k_inner = 0
        while (k_inner ≤ iter_inner) && (f > f₀ + tol_ls * slope * 0.5^k_inner + 1e-6)
            k_inner += 1
            y    .= x₀ + d * 0.5^k_inner         # [n,1]
            X    .= R + diagm(y)                 # [n,n]
            λ, P  = eigen(X)                     # [n,1], [n,n]
            λ, P = Vector{T}(λ), Matrix{T}(P)
            λ    .= reverse(λ)                   # [n,1], [n,n]
            P    .= reverse(P, dims=2)           # [n,n]
            f, Fy = npd_gradient(y, λ, P, b₀, n) # [1], [n,1]
        end

        x₀  = copy(y) # [n,1]
        f₀  = f

        X       .= npd_pca(b₀, X, λ, P, n) # [n,n]
        val_dual = val_R - f₀
        val_obj  = 0.5 * norm(X - R)^2
        gap      = (val_obj - val_dual) / (1 + abs(val_dual) + abs(val_obj))
        b        = b₀ - Fy
        norm_b   = norm(b)
        norm_b_rel      = norm_b / norm_b0

        Ω₀ = npd_set_omega(λ, n) # [n,n] or [r,s]

        k += 1
    end

    X[diagind(X)] .+= τ
    return cov2cor(X)
end

#=
    Return f(yₖ) and ∇f(yₖ) where

    ```math
    f(y) = \\frac{1}{2} \\Vert (A + diag(y))_+ \\Vert_{F}^{2} - e^{T}y
    ```

    and

    ```math
    \\nabla f(y) = Diag((A + diag(y))_+) - e
    ```
=#
function npd_gradient(y::Vector{T}, λ₀::Vector{T}, P::Matrix{T}, b₀::Vector{T}, n::Int) where {T<:AbstractFloat}
    r = sum(λ₀ .> 0)
    λ = copy(λ₀)

    if r == 0
        return (zero(T), zeros(T, n))
    else
        λ[λ .< 0] .= zero(T)
        Fy = Vector{T}(vec(sum((P .* λ') .* P, dims=2)))
        f  = T(T(0.5) * sum(λ.^2) - sum(b₀ .* y))
        return (f, Fy)
    end
end

function npd_pca(b::Vector{T}, X::Matrix{T}, λ::Vector{T}, P::Matrix{T}, n::Int) where {T<:AbstractFloat}
    r = sum(λ .> 0)
    s = n - r

    if r == 0
        X .= zeros(T, n, n)
    elseif r == n
        nothing
    elseif r == 1
        X .= (λ[1] * λ[1]) * (P[:,1] * P[:,1]')
    elseif r ≤ s
        P₁   = @view P[:, 1:r]
        λ₁   = sqrt.(λ[1:r])
        P₁λ₁ = P₁ .* λ₁' # each row of P₁ times λ₁
        X .= P₁λ₁ * P₁λ₁'
    else
        P₂   = @view P[:, (r+1):n]
        λ₂   = sqrt.(-λ[(r+1):n])
        P₂λ₂ = P₂ .* λ₂' # each row of P₂ times λ₂
        X .= X .+ P₂λ₂ * P₂λ₂'
    end

    d  = diag(X)
    d .= max.(d, b)
    X[diagind(X)] .= d
    d .= sqrt.(b ./ d)
    d₂ = d * d'
    X .= X .* d₂
    X
end

# Preconditioned conjugate gradient method to solve Vₖdₖ = -∇f(yₖ)
function npd_pre_cg(
    b::Vector{T},
    c::Vector{T},
    Ω₀::Matrix{T},
    P::Matrix{T},
    ϵ::Real,
    N::Int,
    n::Int) where {T<:AbstractFloat}

    ϵ_b = T(ϵ) * norm(b)

    r   = copy(b)
    z   = r ./ c
    d   = copy(z)
    rz1 = sum(r .* z)
    rz2 = one(T)
    p   = zeros(T, n)
    w   = zeros(T, n)

    for k in 1:N
        if k > 1
            d .= z + d * (rz1 / rz2)
        end

        w .= npd_jacobian(d, Ω₀, P, n)

        denom = sum(d .* w)
        normr = norm(r)

        denom ≤ 0 && return Vector{T}(d / norm(d))

        α = rz1 / denom
        p .+= α*d
        r .-= α*w

        norm(r) ≤ ϵ_b && return Vector{T}(p)

        z .= r ./ c
        rz2, rz1 = copy(rz1), sum(r .* z)
    end

    Vector{T}(p)
end

# Create the precondition matrix used in solving the linear system
# Vₖdₖ = -∇f(yₖ) in the conjugate gradient method.
function npd_precond_matrix(Ω₀::Matrix{T}, P::Matrix{T}, n::Int) where {T<:AbstractFloat}
    r, s = size(Ω₀)

    r == 0 || r == n && return ones(T, n)

    H  = adjoint(P .* P)
    H₁ = @view H[1:r,:]
    H₂ = @view H[r+1:n,:]

    if r < s
        H12 = H₁' * Ω₀
        c   = sum(H₁, dims=1)'.^2 + 2 * sum(H12 .* H₂', dims=2)
    else
        H12 = (1.0 .- Ω₀) * H₂
        c   = sum(H, dims=1)'.^2 - sum(H₂, dims=1)'.^2 - 2 * sum(H₁ .* H12, dims=1)'
    end

    c[c .< 1e-8] .= 1e-8
    Vector{T}(vec(c))
end

# Used in creating the precondition matrix.
function npd_set_omega(λ::Vector{T}, n::Int) where {T<:AbstractFloat}
    r = sum(λ .> 0)
    s = n - r

    r == 0 && return zeros(T, 0, 0)
    r == n && return ones(T, n, n)

    M = zeros(T, r, s)
    λᵣ = @view λ[1:r]
    λₛ = @view λ[r+1:n]
    @inbounds for j = 1:s, i = 1:r
        M[i,j] = λᵣ[i] / (λᵣ[i] - λₛ[j])
    end

    Matrix{T}(M)
end

# Create the Generalized Jacobian matrix for the Newton direction step.
function npd_jacobian(x::Vector{T}, Ω₀::Matrix{T}, P::Matrix{T}, n::Int) where {T<:AbstractFloat}

    r, s = size(Ω₀)
    perturbation = 1e-10

    r == 0 && return zeros(T, n)
    r == n && return Vector{T}(x .* (1 + perturbation))

    P₁ = @view P[:, 1:r]
    P₂ = @view P[:, r+1:n]

    if r < s
        H₁ = diagm(x) * P₁
        Ω  = Ω₀ .* (H₁' * P₂)

        HT₁ = P₁ * P₁' * H₁ + P₂ * Ω'
        HT₂ = P₁ * Ω

        return Vector{T}(vec(sum(P .* [HT₁ HT₂], dims=2) + x .* perturbation))
    else
        H₂ = diagm(x) * P₂
        Ω  = (1 .- Ω₀) .* (P₁' * H₂)

        HT₁ = P₂ * Ω'
        HT₂ = P₂ * H₂' * P₂ + P₁ * Ω

        return Vector{T}(vec(x .* (1 + perturbation) - sum(P .* [HT₁ HT₂], dims=2)))
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -3.0943792528344187 <= q <= 2.261396294653902)

                     add_NL_constraint(m, :(log(1 + exp(-0.6000573802834048 + -0.4460703798694121*$(x[1]) + -0.8919582737388767*$(x[2]))) + log(1 + exp(0.4850728208540551 + 0.8239401438043505*$(x[1]) + -0.080858577965925*$(x[2]))) + log(1 + exp(-0.3015069196609086 + 0.14649274768979348*$(x[1]) + 0.2885676506758026*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module SeisPlot

using PyCall, Statistics, Requires, FFTW
using PyPlot
import SeisMain

#@require PyPlot
const pl = PyNULL()
const lines= PyNULL()

export SeisPlotTX,
SeisPlotAmplitude,
SeisPlotCoordinates,
SeisPlotFK

include("SeisPlotTX.jl")
include("SeisPlotCoordinates.jl")
include("SeisPlotAmplitude.jl")
include("SeisPlotFK.jl")

export pl, lines

function __init__()
    copy!(pl,pyimport_conda("matplotlib.pyplot","matplotlib"))
    copy!(lines,pyimport_conda("matplotlib.lines","matplotlib"))
end
end