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

using Markdown
using InteractiveUtils

# ╔═╡ 6a66d77c-17a4-11eb-0061-f3703df5864c
begin
	using DataFrames
	using CSV
	using Plots
	using GLM
	using Statistics
	using StatsModels
	using StatsPlots
	using Latexify
	using Distributions
end

# ╔═╡ c98883c0-1794-11eb-2d9b-4138e0d901d8
md"""_Test 5_

Submission by: _**Zhuofei, Zhou**_(fregot@icloud.com)

_julia version: 1.5.2_
"""

# ╔═╡ fab7c734-1793-11eb-2c22-3f7a797bd617
md"# Modification of model inapplicability transformation and weighting"

# ╔═╡ 30b3cf96-17a0-11eb-3e5b-cd68a544c808
md"here,import some packages:"

# ╔═╡ f71e2f66-17ac-11eb-18cc-e1ae7832fecc
md"some functions:"

# ╔═╡ ff14ac7a-17ac-11eb-2385-c3565c4ed946
function calculate_one_R2(x::Vector,y::Vector)
    m = size(x)[1]
    X = Float64[x ones(m)]
    β = inv(X'*X)*X'*y
    SSc = y - X*β
    dfc = m - 2
    MSc = SSc/dfc
    e = SSc
    SSh = X*β - mean(y)*ones(m)
    SSH = SSh'*SSh
    SSC = SSc'*SSc
    SST = SSC + SSH
    R = SSH/SST*100
    return R
end

# ╔═╡ 08c6a4ba-17af-11eb-35ff-2f7daa814291

function normal_one_probility(x, y)
    m = size(x)[1]
    X = Float64[x ones(m)]
    β = inv(X'*X)*X'*y
    SSc = y - X*β
    SSC = SSc'*SSc
    dfc = m - 2
    MSc = SSC/dfc
    e = X*β - y
    H = X*inv(X'*X)*X'
    Sᵢ = zeros(m)
    tᵢ = zeros(m)
    p = zeros(m)
    rᵢ = zeros(m)
    press = zeros(m)
    for i  in 1:m
        Sᵢ[i] = ((m-1)*MSc - e[i]/(1 - H[i,i]))/(m - 2)
        tᵢ[i] = e[i]/sqrt(Sᵢ[i]^2*(1 - H[i, i]))
        rᵢ[i] = e[i]/(MSc*(1 - H[i, i]))
        p[i] = (i - 0.5)/m
        press[i] = (e[i]/(1-H[i, i]))^2
    end
    res = DataFrame(y_pre = X*β, e = SSc, Sᵢ = Sᵢ,rᵢ = rᵢ, Tᵢ = tᵢ, p = p, PRESS = press)
    #println(res)
    return res
end

# ╔═╡ ec6c2cc6-17b4-11eb-0b5e-a90d68d8c6b0
function NormalProbabilityPlot(data)
	mu=mean(data)
	sig=std(data)
	n=length(data)
	p=[(i -0.5)/n for i in 1:n]
	x=quantile(Normal(), p)
	y=sort([(i-mu)/sig for i in data])
	qqplot(x, y, qqline=:R, xlabel="Theoretical Quantiles", ylabel="Standardized Residuals", title="Normal Q-Q")
end

# ╔═╡ b8ce3ebc-18f0-11eb-3dd5-3f799465f5ec
correct(text) = Markdown.MD(Markdown.Admonition("correct", "Ans", [text]))

# ╔═╡ 3ef4c8bc-18f4-11eb-3f4a-7179d5e3f071
note(text) = Markdown.MD(Markdown.Admonition("note", "Note", [text]))

# ╔═╡ c7bd8dae-1794-11eb-145e-4938a72ebf82
md"""

## 5.1
Byers and Williams investigated the effect of temperture(regression variable) on the viscosity (response variable) of toluene-tetrahydronaphthalene mixtures.
> Byers and williams研究了温度（回归变量）对甲苯-四氢化萘混合物黏度（响应变量）的影响
1. make a scatter diagram to see if the linear relationship is reasonable.
2. fit the linear model, calculate the President messurement, and make the residual graph. What conclusions can be drawn about the applicable performance of the model
3. The basic priciples of physical chemistry indicate that viscosity is an exponential function of tempterature, redo 2 using transformation.

Data:
"""

# ╔═╡ 20e05c08-17a0-11eb-17d2-c3862d4f222c
begin
	t = [24.9, 35.0, 44.9, 55.1, 65.2, 75.2, 85.2, 95.2]
	y = [1.133, 0.9722, 0.8532, 0.7550, 0.6723, 0.6021, 0.5420, 0.5074]
	data = DataFrame(t=t, y=y)
end

# ╔═╡ 8e0208d6-17a5-11eb-3571-8b4708ecc3de
md"Solutions:
"

# ╔═╡ a9b16f9a-17a5-11eb-1c8b-7b6b4e5bdc68
scatter(t, y, xlabel="t/temperture", ylabel="Pa⋅S", leg=:none)

# ╔═╡ 2d8d245a-17ac-11eb-0d4c-dfb10eeb06e6
ols1 = lm(@formula(y ~ t), data)

# ╔═╡ 26ebed02-17ae-11eb-37f3-8b945ab3fed9
begin
	e = y - predict(ols1)
	scatter(predict(ols1),e)
end

# ╔═╡ 6fff184e-17ad-11eb-3801-0b72654cb699
md"""
$ R^2 = $
"""

# ╔═╡ 4e4dc68c-17ad-11eb-09f1-4fad10ca45e4
latexify(calculate_one_R2(t,y); fmt = "%.2f")

# ╔═╡ 8dbe1f26-17ae-11eb-0f57-c77e8253bb96
correct(md"the residual plot shows a nonlinear pattern and normality is violated.")

# ╔═╡ 908a0072-17af-11eb-07fa-9938001bde9b
res = normal_one_probility(t, y)

# ╔═╡ f650b660-17b1-11eb-065c-f35472d4cf72
res.Tᵢ

# ╔═╡ bae36286-17b0-11eb-0740-53d8bb54e0cc
scatter(t,[log.(y),y])

# ╔═╡ d53f976a-17b0-11eb-147b-c193db7df691
res1 = normal_one_probility(t, log.(y))

# ╔═╡ 03f6bdf4-17b1-11eb-15bb-6d1d0ae5be04
scatter(res1.y_pre, res1.rᵢ)

# ╔═╡ 9f737ec4-17b4-11eb-1384-351e209c6cb8
NormalProbabilityPlot(res.e)

# ╔═╡ a8b71038-17b2-11eb-3ac4-d5ff9fffaa11
correct(md"There is slight improvement in the model")

# ╔═╡ fe82916a-18f0-11eb-0f29-2fa6d9326c13
html"

"

# ╔═╡ 095bae32-18f1-11eb-3ff7-cfdc4125d745
md"""
## 5.2

The vapor pressure of water at different temperatures.

1. make a scatter diagram to see if the linear relationship is reasonable.

2. fitting linear model, and calculate President messurement, and make the residual graph. What conclusions can be drawn about the applicable performance of the model.

3.  The clausius-Clapeyron equation in physical chemistry states: $ln(p_c) \sim -\frac1T$, based on the information, do some appropriate transformission.
"""

# ╔═╡ 8f80cc44-18f2-11eb-12d4-5727ca38789d
#data
begin
	mmHg = [4.6, 9.2, 17.5, 31.8, 55.3, 92.5, 149.4, 233.7, 355.1, 525.8, 760.0]
	data2 = DataFrame(t = [273, 283, 293, 303, 313, 323, 333, 343, 353, 363, 373], mmHg = mmHg)
end

# ╔═╡ 029f2540-18f3-11eb-3523-f35f3e2ac61f
scatter(data2.t, data2.mmHg, leg=:none, xlabel="temperatures", ylabel="mmHg")

# ╔═╡ 383465bc-18f3-11eb-09b6-51709b7a4053
correct(md"absolutely, linear relationship is unreasonable.")

# ╔═╡ 05a7a704-18f4-11eb-3bf0-3de804963e50
ols2= lm(@formula(mmHg ~ t), data2)

# ╔═╡ 595e4ea4-18f3-11eb-3e10-8f972926492d
begin
	res2 = normal_one_probility(data2.t, data2.mmHg)
	scatter(res2.y_pre, res2.e, title="residuals plot", leg=:none)
end

# ╔═╡ c14304bc-18f3-11eb-087a-4fe3c1c9f1ce
correct(md"Direct use of linear model isn't a good way.")

# ╔═╡ f27b2e74-18f3-11eb-1978-73bd712aec51
note(md"we are kown ``ln(p) \sim -\frac1T``,so we can do like:``y' = ln(mmHg), x'=\frac1T``")

# ╔═╡ a5642522-18f4-11eb-1c36-b3e1e39249b1
begin
	ŷ₂ = log.(data2.mmHg)
	x̂₂ = 1 ./ data2.t
	data2.y_t = ŷ₂
	data2.x_t = x̂₂
	ols22 = lm(@formula(y_t ~ x_t), data2)
end

# ╔═╡ 5f58147a-18f5-11eb-2004-75b188ce8a86
begin
	res22 = normal_one_probility(data2.x_t, data2.y_t)
	scatter(res22.y_pre, res22.e, title="after transmission residuals plot", leg=:none)
end

# ╔═╡ 89c38532-18f5-11eb-2b17-0f367a4a51fd
scatter(data2.x_t, data2.y_t, leg=:none)

# ╔═╡ ae63ac1e-18f5-11eb-2790-c5d647e2c52a
correct(md"the function: ``ln(mmHg) = -5200.76 \cdot \frac1T + 20.6074``")

# ╔═╡ dbe4f300-18f5-11eb-3751-2f7a7997a159
html"

"

# ╔═╡ b5499d0a-1929-11eb-346c-2fc9cda07625
md"""
## 5.3
The average number of living 🦠(bacteria) in canned food products and the number of minutes exposed to 300∘F.

1. make a scatter diagram to see if the linear relationsihp is reasionable.

2. fitting linear model, calculate President messurement, and make the residual graph. What conclusions can be drawn about the applicable performance of the model.

3. The suitable model transformation for this data is identified, and the data is fitted with the transformed model and the general test of model applicability is carried out.
"""

# ╔═╡ a85eb29c-192b-11eb-0d05-c1087fb1ee63
data3 = DataFrame(time=[1,2,3,4,5,6,7,8,9,10,11,12],number=[175,108,95,82,71,50,49,31,28,17,16,11])

# ╔═╡ f5672d3a-192b-11eb-271d-3f82c74c1fce
scatter(data3.time, data3.number, leg=:none, xlabel="expose time", ylabel="bacteria numbers")

# ╔═╡ 80c4ddf0-192c-11eb-02eb-6b8e5ed1ecde
correct(md"Direct use of linear model isn't a good way.")

# ╔═╡ 88253dec-192c-11eb-3fe4-d3c9f3eb9882
ols3 = lm(@formula(number ~ time), data3)

# ╔═╡ f1750ab0-192d-11eb-23c4-f9f82e4343e7
r2(ols3)

# ╔═╡ a94f5610-192c-11eb-387f-87f04e2d3774
begin
	res3 = normal_one_probility(data3.time, data3.number)
	scatter(res3.y_pre, res3.e, leg=:none, xlabel="predict y value", ylabel="residual")
end

# ╔═╡ b29e4abc-192e-11eb-23ce-c52e1b6c41e0
NormalProbabilityPlot(res3.e)

# ╔═╡ c0f3ad04-192f-11eb-166b-59f0ccd50ce1
note(md"we can find a obvious outliers point, if we delect the point.

here we choose two ways to do linear transformation.

- ``y=\beta_0 e ^{\beta_1 x}``
- ``y = \beta_0 + \beta_1 log x`` ")

# ╔═╡ 9ac087a0-1930-11eb-28e9-e7886e5816bc
scatter(data3.time, log.(data3.number))

# ╔═╡ b367c318-1930-11eb-0d0b-53498dd6eab0
scatter(log.(data3.time), data3.number)

# ╔═╡ aceb07f2-1933-11eb-122b-4185d4d50031
begin
	res33 = normal_one_probility(log.(data3[2:end, :].time), data3[2:end, :].number)
	scatter(res33.y_pre, res33.e, leg=:none)
end

# ╔═╡ 00a65f8a-1941-11eb-0ece-7beb73a126d6
correct(md"so, if we delect a outlier point, and do a ``x' = ln(x)``")

# ╔═╡ a74d3a02-1941-11eb-2192-510011d1f284
html"

"

# ╔═╡ dee284c2-1941-11eb-0dab-ab42595e7d56
md"""
## 5.4
Considering the data shown below, the scatter diagram is constructed and the appropriate form of the refression model is suggested.
"""

# ╔═╡ 0102859c-1943-11eb-332e-3f896f4114e5
data4 = DataFrame(x=[10,15,18,12,9,8,11,6], y=[0.17,0.13,0.09,0.15,0.20,0.21,0.18,0.24])

# ╔═╡ 32ab63f2-1943-11eb-2ece-03627d4a4e09
scatter(data4.x, data4.y, leg=:none)

# ╔═╡ 4c06b20c-1943-11eb-3191-21e1fb66dbc8
ols4 = lm(@formula(y ~ x), data4)

# ╔═╡ 69ff0ebc-1943-11eb-101b-43f5f6d7ce18
begin
	res4 = normal_one_probility(data4.x, data4.y)
	scatter(res4.y_pre, res4.e, xlabel="fitting value", ylabel="residual")
end

# ╔═╡ 9afdcdfa-1943-11eb-26e4-eb1268dfe551
correct(md"``y = -0.0121215 \cdot x + 0.306103``, but the residual plot does not look good. Taking the _natural log of x_ makes better model.")

# ╔═╡ 6e64f52e-1944-11eb-05db-49ed550f6fa3
begin
	res44 = normal_one_probility(log.(data4.x), data4.y)
	scatter(res44.y_pre, res44.e, xlabel="fitting value", ylabel="residual", label="do log")
end

# ╔═╡ e3cfa076-1943-11eb-35ac-3110da30870f
html"

"

# ╔═╡ 9d324cb2-1944-11eb-38b1-ff70b0c41486
md"""
## 5.5
Bottle makers have recoverded the average number of defects due to gravel per 10,000 bottles and the number of weeks since the last furnace was repaired.

1. fitting the linear model and Standard test for model suitability
2. suggeat a reasionable transformation to slove the mistake in 1, fitting the after transformation and model suitability
"""

# ╔═╡ c769f1e0-1946-11eb-397e-99b41d51b867
data5 = DataFrame(number=[13.0,16.1,14.5,17.8,22.0,27.4,16.8,34.2,65.6,49.2,66.2,81.2,87.4,114.5], weeks=[4,5,6,7,8,9,10,11,12,13,14,15,16,17])

# ╔═╡ 2215d92e-1947-11eb-1840-4d92bcb14185
scatter(data5.weeks, data5.number, leg=:none, xlabel="weeks", ylabel="number")

# ╔═╡ 14111252-19a2-11eb-16c2-fbb07b4d574f
ols5 = lm(@formula(number ~ weeks), data5)

# ╔═╡ 2e46c27a-19a2-11eb-31f6-cb94419d9a62
begin
	res5 = normal_one_probility(data5.weeks, data5.number)
	scatter(res5.y_pre, res5.e, xlabel="fitting value", ylabel="residual")
end

# ╔═╡ 54c592f0-19a2-11eb-0126-21695f895229
NormalProbabilityPlot(res5.e)

# ╔═╡ 87fe612e-19a2-11eb-1513-a3aafbf41113
begin
	res55 = normal_one_probility(data5.weeks, log.(data5.number))
	plot(scatter(res55.y_pre, res55.e, xlabel="fitting value", ylabel="residual", leg=:none), NormalProbabilityPlot(res55.e))
end

# ╔═╡ 46694f30-1b91-11eb-269f-6dede05ef8fb
scatter(data5.weeks, log.(data5.number), leg=:none, xlabel="weeks", ylabel="log(numbers)")

# ╔═╡ f87c19ba-1b92-11eb-295c-41e43149192c
begin
	change_data5 = DataFrame(weeks = data5.weeks, numbers = log.(data5.number))
	ols55 = lm(@formula(numbers ~ weeks), change_data5)
end

# ╔═╡ 47b98d3a-1b93-11eb-12c1-1d21ee1f2305
correct(md"""so the function: $ln(numbers) = 1.71622 + 0.1735\cdot weeks$""")

# ╔═╡ 839c68d6-1b93-11eb-3fef-0db10bfe7575
html"

"

# ╔═╡ d389d55e-1b93-11eb-0664-e1438da2e017
md"""
## 5.6
Consider the fuel consumpsion data in table B-18, ignoring the regreesion variable x₁, and recalling the complete residual analysis done in exercise 4.27, can transformation improve this analysis? why or why not?

> 考虑表B-18中的燃料消耗数据，忽略回归变量x1，回忆习题4.27所做的完整的残差分析，变换能改善这一分析吗？为什么能或为什么不能？
"""

# ╔═╡ 28e28770-1b95-11eb-0f93-b7709ed092d4
html"

"

# ╔═╡ 30098fc6-1b95-11eb-362e-bddf13d1e65b
md"""
## 5.7
Consider the methanol oxidation data in Table B-20 and make a throunth analysis of this data.Recalling the residual analysis of this data in problem 4.29, Can transformation improve this analysis? why or why not?

>考虑表B-20中的甲醇氧化数据，对这一数据做彻底的分析，回忆习题4.29对这一数据所做的残差分析。变换能改善这一分析吗？为什么能或为什么不能？

"""

# ╔═╡ b8825eaa-1b95-11eb-2667-3d5b149e2792
html"

"

# ╔═╡ c4edfb5e-1b95-11eb-06cc-a3deca511348
md"""
## 5.8

Consider three models

1. $y = \beta_0 + \beta_1 ( 1/x) + \epsilon$
2. $1/y = \beta_0 + \beta_1x + \epsilon$
3. $y = x/(\beta_0 - \beta_1 x) + \epsilon$

All three of these models are linearized by reciprocal transformation, sketch out what features appear in the scatter plot that will allow you to pick one of the three models.
> 这三个模型都通过倒数变换而线性化， 画出草图， 散点图出现什么特征会让你挑选三种模型中的一种。
"""

# ╔═╡ bd9932e6-1b96-11eb-35c7-957ea3fef253
html"

"

# ╔═╡ Cell order:
# ╟─c98883c0-1794-11eb-2d9b-4138e0d901d8
# ╟─fab7c734-1793-11eb-2c22-3f7a797bd617
# ╟─30b3cf96-17a0-11eb-3e5b-cd68a544c808
# ╟─6a66d77c-17a4-11eb-0061-f3703df5864c
# ╟─f71e2f66-17ac-11eb-18cc-e1ae7832fecc
# ╟─ff14ac7a-17ac-11eb-2385-c3565c4ed946
# ╟─08c6a4ba-17af-11eb-35ff-2f7daa814291
# ╟─ec6c2cc6-17b4-11eb-0b5e-a90d68d8c6b0
# ╟─b8ce3ebc-18f0-11eb-3dd5-3f799465f5ec
# ╟─3ef4c8bc-18f4-11eb-3f4a-7179d5e3f071
# ╟─c7bd8dae-1794-11eb-145e-4938a72ebf82
# ╟─20e05c08-17a0-11eb-17d2-c3862d4f222c
# ╟─8e0208d6-17a5-11eb-3571-8b4708ecc3de
# ╟─a9b16f9a-17a5-11eb-1c8b-7b6b4e5bdc68
# ╟─2d8d245a-17ac-11eb-0d4c-dfb10eeb06e6
# ╟─26ebed02-17ae-11eb-37f3-8b945ab3fed9
# ╟─6fff184e-17ad-11eb-3801-0b72654cb699
# ╟─4e4dc68c-17ad-11eb-09f1-4fad10ca45e4
# ╟─8dbe1f26-17ae-11eb-0f57-c77e8253bb96
# ╠═908a0072-17af-11eb-07fa-9938001bde9b
# ╠═f650b660-17b1-11eb-065c-f35472d4cf72
# ╠═bae36286-17b0-11eb-0740-53d8bb54e0cc
# ╠═d53f976a-17b0-11eb-147b-c193db7df691
# ╠═03f6bdf4-17b1-11eb-15bb-6d1d0ae5be04
# ╟─9f737ec4-17b4-11eb-1384-351e209c6cb8
# ╟─a8b71038-17b2-11eb-3ac4-d5ff9fffaa11
# ╟─fe82916a-18f0-11eb-0f29-2fa6d9326c13
# ╟─095bae32-18f1-11eb-3ff7-cfdc4125d745
# ╟─8f80cc44-18f2-11eb-12d4-5727ca38789d
# ╟─029f2540-18f3-11eb-3523-f35f3e2ac61f
# ╟─383465bc-18f3-11eb-09b6-51709b7a4053
# ╠═05a7a704-18f4-11eb-3bf0-3de804963e50
# ╠═595e4ea4-18f3-11eb-3e10-8f972926492d
# ╟─c14304bc-18f3-11eb-087a-4fe3c1c9f1ce
# ╟─f27b2e74-18f3-11eb-1978-73bd712aec51
# ╟─a5642522-18f4-11eb-1c36-b3e1e39249b1
# ╟─5f58147a-18f5-11eb-2004-75b188ce8a86
# ╟─89c38532-18f5-11eb-2b17-0f367a4a51fd
# ╟─ae63ac1e-18f5-11eb-2790-c5d647e2c52a
# ╟─dbe4f300-18f5-11eb-3751-2f7a7997a159
# ╟─b5499d0a-1929-11eb-346c-2fc9cda07625
# ╟─a85eb29c-192b-11eb-0d05-c1087fb1ee63
# ╠═f5672d3a-192b-11eb-271d-3f82c74c1fce
# ╟─80c4ddf0-192c-11eb-02eb-6b8e5ed1ecde
# ╠═88253dec-192c-11eb-3fe4-d3c9f3eb9882
# ╠═f1750ab0-192d-11eb-23c4-f9f82e4343e7
# ╠═a94f5610-192c-11eb-387f-87f04e2d3774
# ╠═b29e4abc-192e-11eb-23ce-c52e1b6c41e0
# ╟─c0f3ad04-192f-11eb-166b-59f0ccd50ce1
# ╠═9ac087a0-1930-11eb-28e9-e7886e5816bc
# ╠═b367c318-1930-11eb-0d0b-53498dd6eab0
# ╠═aceb07f2-1933-11eb-122b-4185d4d50031
# ╟─00a65f8a-1941-11eb-0ece-7beb73a126d6
# ╟─a74d3a02-1941-11eb-2192-510011d1f284
# ╟─dee284c2-1941-11eb-0dab-ab42595e7d56
# ╟─0102859c-1943-11eb-332e-3f896f4114e5
# ╟─32ab63f2-1943-11eb-2ece-03627d4a4e09
# ╟─4c06b20c-1943-11eb-3191-21e1fb66dbc8
# ╟─69ff0ebc-1943-11eb-101b-43f5f6d7ce18
# ╟─9afdcdfa-1943-11eb-26e4-eb1268dfe551
# ╟─6e64f52e-1944-11eb-05db-49ed550f6fa3
# ╟─e3cfa076-1943-11eb-35ac-3110da30870f
# ╟─9d324cb2-1944-11eb-38b1-ff70b0c41486
# ╟─c769f1e0-1946-11eb-397e-99b41d51b867
# ╟─2215d92e-1947-11eb-1840-4d92bcb14185
# ╟─14111252-19a2-11eb-16c2-fbb07b4d574f
# ╟─2e46c27a-19a2-11eb-31f6-cb94419d9a62
# ╟─54c592f0-19a2-11eb-0126-21695f895229
# ╟─87fe612e-19a2-11eb-1513-a3aafbf41113
# ╟─46694f30-1b91-11eb-269f-6dede05ef8fb
# ╟─f87c19ba-1b92-11eb-295c-41e43149192c
# ╟─47b98d3a-1b93-11eb-12c1-1d21ee1f2305
# ╟─839c68d6-1b93-11eb-3fef-0db10bfe7575
# ╟─d389d55e-1b93-11eb-0664-e1438da2e017
# ╟─28e28770-1b95-11eb-0f93-b7709ed092d4
# ╟─30098fc6-1b95-11eb-362e-bddf13d1e65b
# ╟─b8825eaa-1b95-11eb-2667-3d5b149e2792
# ╟─c4edfb5e-1b95-11eb-06cc-a3deca511348
# ╟─bd9932e6-1b96-11eb-35c7-957ea3fef253

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_name: ZhouZhuofei/juliaDL


text:
function read_credit_spread_data(dt::Date)

	url = "https://www.anbima.com.br/informacoes/curvas-debentures/CD-down.asp"
	headers = [
	    "Host" => "www.anbima.com.br"
	    "User-Agent" => "Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:80.0) Gecko/20100101 Firefox/80.0"
	    "Accept" => "text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,*/*;q=0.8"
	    "Accept-Language" => "pt-BR,pt;q=0.8,en-US;q=0.5,en;q=0.3"
	    "Accept-Encoding" => "gzip, deflate, br"
	    "Content-Type" => "application/x-www-form-urlencoded"
	    "Origin" => "https://www.anbima.com.br"
	    "Connection" => "keep-alive"
	    "Referer" => "https://www.anbima.com.br/informacoes/curvas-debentures/default.asp"
	    "Upgrade-Insecure-Requests" => "1"
	]
	str_dt = Dates.format(dt, "dd%2Fmm%2Fyyyy")
	body = "escolha=2&Idioma=PT&saida=csv&Dt_Ref_Ver=20200918&Dt_Ref=$str_dt"

	r = HTTP.post(url, headers, body)
	response_date = _parse_date_credit_curve(r)
	@assert dt == response_date "Wrong date: expected $dt, got $response_date"

	response_body = String(r.body)
	@assert _has_data_credit_curve(response_body) "no credit spread data for date $dt"
	out = _parse_credit_curve_result(response_body)
	return out

end

function _has_data_credit_curve(str::String)
	spl = split(str, "\r\n")
	return length(spl) > 2
end

function _parse_date_credit_curve(r)
	d = Dict(r.headers)
	@assert haskey(d, "content-disposition") "response header must have content-disposition"
	m = match(r"\d{8}", d["content-disposition"])
	@assert m != nothing "invalid content-disposition"
	str_dt = String(m.match)

	return Date(str_dt, "ddmmyyyy")
end

function _parse_credit_curve_result(str::String)
	rates = _parse_credit_curve_rates(str)
	return rates
end

function _parse_credit_curve_rates(str::String)
	lines = split(str, "\r\n")
	vertices = Vector{Float64}()
	rate_aaa = Vector{Float64}()
	rate_aa  = Vector{Float64}()
	rate_a   = Vector{Float64}()
	for line in lines[3:end]
		if line != ""
			s = split(line, ";")
			push!(vertices, parse_value(Float64, s[1]))
			push!(rate_aaa, parse_value(Float64, s[2]))
			push!(rate_aa , parse_value(Float64, s[3]))
			push!(rate_a  , parse_value(Float64, s[4]))
		end
	end
	return Dict([
		"AAA" => [k => v/100 for (k,v) in zip(vertices, rate_aaa)]
		"AA" =>  [k => v/100 for (k,v) in zip(vertices, rate_aa)]
		"A" =>   [k => v/100 for (k,v) in zip(vertices, rate_a)]
	])
end

function calibrate_credit_spread_curves(dt::Date)
	data = read_credit_spread_data(dt)
	return _do_calibrate(data)
end

function _do_calibrate(d::Dict{String, Vector{Pair{Float64, Float64}}}; best_of=3)

	taus = Dict{String, Vector{Int64}}()
	rates = Dict{String, Vector{Float64}}()

	for rating in keys(d)
		taus[rating] = [Int64(round(252*first(x))) for x in d[rating]]
		rates[rating] = [last(x) for x in d[rating]]
	end

	function get_curves(x)
		@assert length(x) == 5
		β1_AAA    = x[1]
		β1_AA     = x[2]
		β1_A      = x[3]
		β_credito = x[4]
		β3 = 0.0  # no curvature
		λ_credito = x[5]
		cvs = Dict([
			"AAA" => NelsonSiegelCurve(β1_AAA, β_credito, β3, λ_credito)
			"AA"  => NelsonSiegelCurve(β1_AA , β_credito, β3, λ_credito)
			"A"   => NelsonSiegelCurve(β1_A  , β_credito, β3, λ_credito)
		])
		return cvs
	end

	function f_optim(x)
		cvs = get_curves(x)
		out = 0.0
		for rating in keys(cvs)
			model_rates = [zerorate(cvs[rating], t) for t in taus[rating]]
			squared_diff = sum((100*model_rates - 100*rates[rating]).^2)
			out += squared_diff
		end
		return out
	end

	best_minimum = Inf
	solution = nothing
	for i in 1:best_of
		#@info "Solution $i/$best_of"
		this_solution = Evolutionary.optimize(
			f_optim,
			rand(5),
			Evolutionary.GA(
				populationSize = 1000,
				crossoverRate=0.5,
				mutationRate=0.5,
				ɛ=0.1,
				selection = susinv,
				crossover = discrete,
				mutation = domainrange(ones(5))
			),
			Evolutionary.Options(
				show_trace=false,
				show_every=100)
		)
		if Evolutionary.converged(this_solution)
			this_min = Evolutionary.minimum(this_solution)
			if this_min < best_minimum
				#@info "BEST!"
				best_minimum = this_min
				solution = this_solution
			end
		end
	end

	if solution == nothing
		error("could not solve spread curve minimization")
	end
	x_optim = Evolutionary.minimizer(solution)
	curve_optim = get_curves(x_optim)

	return curve_optim

end

ext: jl


text:
mutable struct Tree
    parent
    num_children
    children
    _size
end

function Tree()
    parent = nothing
    num_children = 0
    children = []
    _size = nothing
    return Tree(parent, num_children, children, _size)
end

function String_tree( t::Tree, level=0)
    for child in t.children
        if child isa typeof(t)
            ret = ""
            ret *= String_tree(child, level+1)
        else
            ret *= "\t"*string(level) * string(child) * "\n"
        end
    end
    return ret

end

function add_child(t::Tree,c)
    if c isa typeof(t)
        c.parent = t
    end
    push!(t.children, c)
    t.num_children = t.num_children + 1
end

function sizeof(t::Tree)
    if t._size != nothing
        return t._size
    end
    size =1
    for i in 1:(t.num_children)
        size = size + length(self.children[i])
    end
    t._size = size
    return size
end

function children_vector(t::Tree)
    r_list = []
    for i in 1:(t.num_children)
        if t.children[i] isa typeof(t)

            push!(r_list, 4)
        else
            push!(r_list, self.children[i])
        end
    end
    return r_list
end

function to_string(t::Tree)
    r_list = []
    for i in 1:(t.num_children)
        if t.children[i] isa typeof(t)
            push!(r_list, "( " * to_string(t.children[i]) * " )")
        else
            push!(r_list, string(t.children[i]))
        end
    end
    return join(r_list, " ")
end

function to_list(t::Tree, form_manager)
    r_list = []
    for i in 1:(t.num_children)
        if t.children[i] isa typeof(t)
            push!(r_list, get_symbol_idx(form_manager, "("))
            cl = to_list(t.children[i], form_manager)
            for k in 1:length(cl)
                push!(r_list, cl[k])
            end
            push!(r_list, get_symbol_idx(form_manager, ")"))
        else
            push!(r_list, t.children[i])
        end
    end
    return join(r_list, " ")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 54530d68-371e-11eb-198f-09c7ddd312d0
md" ## _Control Flow_
In this lesson we will learn how to work with control statements. We will first learn how to use conditional blocks like `if ... else` blocks and then we will learn how to perform loops.
"

# ╔═╡ eedcadc6-371e-11eb-157d-519a4216048e
md" ## If...else
When a program needs to take decisions according to certain conditions, the `if ... else` block is the default choice.

Let’s suppose that we want to write a simple implementation of the absolute value of a number. The absolute value of a number is defined as the number itself, if the number is positive, or the opposite if the value is negative. This is the typical case where the `if ... else` construct is useful! We can write a simple absolute function in this way:
"

# ╔═╡ 0d8b49a8-371f-11eb-2095-f396ca04dbf8
function absolute(x)
    if x >= 0
        return x
    else
        return -x
    end
end

# ╔═╡ 141b93c2-371f-11eb-32d6-27f16f613bb5
md" As you can see, an `if ... else` block is closed with the word `end`, like a function.

If we need to check more than one condition, we can bind together conditions using:

- “and” is written as `&` (if both statements are true return true, else return false)
- “or” is written as `||` (if at least one statement is true return true, else return false)
For example, if we want to check whether 3 is both minor than 4 and major than 1 we type:
"

# ╔═╡ 5bd9ab0c-371f-11eb-34c4-3bf036e74aa2
if 1 < 3 & 3 < 4
    md"eureka!"
end

# ╔═╡ 7612ee90-371f-11eb-1bd1-1901a20d17ac
md"
If we want to check if a value satisfies one of several different conditions, we can use the elseif statement: Julia will check if the first condition specified in the if is satisfied, if it is not met it moves on to the first elseif and so on.
"

# ╔═╡ 8fad36a8-371f-11eb-2813-799af7e042d6
x=109

# ╔═╡ 661d3b08-371f-11eb-2de6-2d6e6159b11b

	if x<1
	    md"$x < 1"
	elseif x < 3
	    md"$x < 3"
	elseif x < 100
	    md"$x < 100"
	else
	    md"$x is really big!"
	end

# ╔═╡ fed3f85e-3720-11eb-06ee-478f071c7445
md"I take the occasion to introduce **string interpolation**. With the indication `$x` we tell Julia that it must substitute to `$x` the value of `x`, in this case 42. This is particularly useful when we want to print values, or we want to make custom messages:
"

# ╔═╡ 58766288-3722-11eb-3f92-b96dc7699a55
name1 = "Zoe zebra"

# ╔═╡ 662ed338-3722-11eb-10f6-0913ff2bb07f
name2 = "Peppa"

# ╔═╡ 663b5b76-3722-11eb-0279-5b549e2422b7
md"Welcome $name1, this is $name2 :)"

# ╔═╡ 4274730c-3723-11eb-31b0-15f4fb936db5
md" ## Loops
A loop is the operation of repeating the same set of instructions several times. Loops are useful when we want to compute the value of a function over several points, we need to perform some operations on the elements of an array or we need to print the elements of a list.
"

# ╔═╡ 5c29e0d4-3723-11eb-210d-e16125372096
md" **For Loops**
........................................................................................................................................................................

Sometimes we want to iterate over a list of values and perform some operation on each element.

For example let’s suppose we want to print all the squares of the numbers comprised between 1 and 10, we can do it using a `for` loop:
"

# ╔═╡ afc4775e-3723-11eb-28c9-5d9248c25d30
for i in 1:10
    md"i^2"
end

# ╔═╡ 8bf7f11a-3724-11eb-343c-cde9b80cc292
md"
`i` is the variable which contains the data which gets updated at each new cycle (in this case i holds the numbers from 1 to 10 respectively), while `1:10` is a **range**. A range is an iterable object which does exactly what its name suggests: it specifies the range on which the loop has to be performed (in this case 1 to 10).

It is also possible to use the alternative notation for i = 1:10 which is completely equivalent.

Please notice that it is possible to loop not only over ranges (which can also be specified using the `range` function) but also lists (i.e. arrays, tuples, etc).

For example, let’s suppose we have a list of persons and we want to greet all of them, we can do it with the `for` statement:
"

# ╔═╡ c46343c4-3724-11eb-1bdf-8558d7e2c915
persons = ["Alice", "Bob", "Carla", "Daniel"]

# ╔═╡ ffc3cf7e-3724-11eb-31a6-ff2ab6952113
for person in persons
    md"Hello $person, welcome to techytok!"
end

# ╔═╡ 0f213132-3725-11eb-0023-e37fd3aa9933
md"Here instead of a range, we iterate over the elements of persons (i.e. the names of the persons that I want to greet) and in this case person will hold the name of a single person, which changes at each iteration step.
"

# ╔═╡ 86a93f60-3725-11eb-3e83-db6d45c4918a
md"**Breaks**
........................................................................................................................................................................

In the case we want to forcefully interrupt a for loop we can use the break statement, for example:
"

# ╔═╡ c287b7be-3725-11eb-3c5e-1578203b2779
for i in 1:100
    if i>10
        break
    else
   	    println(i^2)
    end
end

# ╔═╡ d4821eb4-3725-11eb-091b-5f9e7fd555ef
md" Here we check if `i>10`, in that case we break the loop and finish the iteration, else we print `i^2`.
"

# ╔═╡ e67787c4-3725-11eb-3659-111648e1efbb
md"**Continue**
........................................................................................................................................................................

This is the opposite of break, continue will forcefully skip the current iteration and move to the next cycle:
"

# ╔═╡ fa813a82-3725-11eb-31c3-915661199563
for i in 1:30
    if i % 3 == 0
        continue
    else
        println(i)
    end
end

# ╔═╡ 0ceee5e0-3729-11eb-3a5a-f58377c6fb9a
md"**While Loop**
........................................................................................................................................................................

When a loop needs to continue until a certain condition is met, a while loop is the preferable choice:
"

# ╔═╡ 1c7f425c-3729-11eb-086a-3b64ef09e50c
function while_test()
    i=0
    while(i<30)
        println(i)
        i += 1
    end
end

# ╔═╡ 234a80b0-3729-11eb-0f94-3b097c8b0260
md" ## Enumerate
`enumerate` is a function which comes in handy when we need to iterate on an array (or similar) and we need to keep track of the number of iterations we have already performed.

`enumerate` will return an iterator (which is something like an array which can be iterated in for loops). It will produce couples of the form `(i, x[i])`.
"

# ╔═╡ 384174f6-3729-11eb-3dca-e7b0d2e6a3b5
z = ["a","b","c"]

# ╔═╡ 3eef4080-3729-11eb-0445-2d68fe408a17
for couple in enumerate(z)
    println(couple)
end

# ╔═╡ 3efc4ce4-3729-11eb-0401-e91c29e82242
(1, "a")

# ╔═╡ 3f405c54-3729-11eb-2792-05f370aa82b0
(2, "b")

# ╔═╡ 3f4e1e68-3729-11eb-0c64-1bf3e2996f25
(3, "c")

# ╔═╡ 4c84b284-3729-11eb-0899-3fd8fdc2b4bd
md"The same result could have been obtained “manually”:"

# ╔═╡ 5b92232e-3729-11eb-3775-611d7fe365f0
y = ["a","b","c"]

# ╔═╡ 63bf58fa-3729-11eb-01a8-3b19208f9751
enum_array = [(1,"a"), (2,"b"), (3,"c")]

# ╔═╡ 63c96604-3729-11eb-0406-01eece93f912
for i in 1:length(y)
    println(enum_array[i])
end

# ╔═╡ 63dc02a4-3729-11eb-212c-5bb44f5ec5cc
(1, "a")

# ╔═╡ 6405dee2-3729-11eb-1a6d-dff129465d17
(2, "b")

# ╔═╡ 6412f142-3729-11eb-1081-d51e4329b48a
(3, "c")

# ╔═╡ Cell order:
# ╟─54530d68-371e-11eb-198f-09c7ddd312d0
# ╟─eedcadc6-371e-11eb-157d-519a4216048e
# ╠═0d8b49a8-371f-11eb-2095-f396ca04dbf8
# ╟─141b93c2-371f-11eb-32d6-27f16f613bb5
# ╠═5bd9ab0c-371f-11eb-34c4-3bf036e74aa2
# ╟─7612ee90-371f-11eb-1bd1-1901a20d17ac
# ╠═8fad36a8-371f-11eb-2813-799af7e042d6
# ╠═661d3b08-371f-11eb-2de6-2d6e6159b11b
# ╟─fed3f85e-3720-11eb-06ee-478f071c7445
# ╠═58766288-3722-11eb-3f92-b96dc7699a55
# ╠═662ed338-3722-11eb-10f6-0913ff2bb07f
# ╠═663b5b76-3722-11eb-0279-5b549e2422b7
# ╟─4274730c-3723-11eb-31b0-15f4fb936db5
# ╠═5c29e0d4-3723-11eb-210d-e16125372096
# ╠═afc4775e-3723-11eb-28c9-5d9248c25d30
# ╟─8bf7f11a-3724-11eb-343c-cde9b80cc292
# ╠═c46343c4-3724-11eb-1bdf-8558d7e2c915
# ╠═ffc3cf7e-3724-11eb-31a6-ff2ab6952113
# ╟─0f213132-3725-11eb-0023-e37fd3aa9933
# ╠═86a93f60-3725-11eb-3e83-db6d45c4918a
# ╠═c287b7be-3725-11eb-3c5e-1578203b2779
# ╟─d4821eb4-3725-11eb-091b-5f9e7fd555ef
# ╟─e67787c4-3725-11eb-3659-111648e1efbb
# ╠═fa813a82-3725-11eb-31c3-915661199563
# ╟─0ceee5e0-3729-11eb-3a5a-f58377c6fb9a
# ╠═1c7f425c-3729-11eb-086a-3b64ef09e50c
# ╟─234a80b0-3729-11eb-0f94-3b097c8b0260
# ╠═384174f6-3729-11eb-3dca-e7b0d2e6a3b5
# ╠═3eef4080-3729-11eb-0445-2d68fe408a17
# ╠═3efc4ce4-3729-11eb-0401-e91c29e82242
# ╠═3f405c54-3729-11eb-2792-05f370aa82b0
# ╠═3f4e1e68-3729-11eb-0c64-1bf3e2996f25
# ╟─4c84b284-3729-11eb-0899-3fd8fdc2b4bd
# ╠═5b92232e-3729-11eb-3775-611d7fe365f0
# ╠═63bf58fa-3729-11eb-01a8-3b19208f9751
# ╠═63c96604-3729-11eb-0406-01eece93f912
# ╠═63dc02a4-3729-11eb-212c-5bb44f5ec5cc
# ╠═6405dee2-3729-11eb-1a6d-dff129465d17
# ╠═6412f142-3729-11eb-1081-d51e4329b48a


text:
# Date: March 1, 2019
# Author: Bernhard König
# Title: An application of Regression Trees for Reserving

# We have used the code in RCode\GenerateIndividualClaimsData.R
# to generate the data Simulated.Cashflows.csv
# We refer to the work of Wuethrich and Gabrielli for details on the
# simulated claims data, see the following references:
    ### Source: https://people.math.ethz.ch/~wmario/simulation.html
    ### Reference: Individual Claims History Simulation Machine. ###
    ###            SSRN Manuscript ID 3130560.                   ###
    ###  Authors: Andrea Gabrielli, Mario V. Wuthrich    ###

#
@assert splitdir(pwd())[2] == "DecisionTrees.jl" """Error. The working directory should be the DecisionTrees.jl folder.\r\nUse 'cd("C:\\Users\\bernhard.konig\\Documents\\ASync\\home\\Code\\Julia\\DecisionTrees.jl")' or similar to change the working dir in Julia """
# if the working dir is not correct, the include command further down will fail

@warn("You need to install DecisionTrees.jl. Consider the readme of this page: 'https://github.com/kafisatz/DecisionTrees.jl'")
@warn("You need to ensure that the required packages are are installed for this code to run.")

using Revise
using Statistics
using StatsBase
import StatsPlots
# using Plots

# these are std lib packages and installed by default
using Distributed, Random, DelimitedFiles,Pkg

Distributed.@everywhere import CSV
Distributed.@everywhere import DataFrames
Distributed.@everywhere import DataFrames: DataFrame,groupby,combine,names!,aggregate

@time Distributed.@everywhere using DecisionTrees # first time usage (precompilation) may take some time

folderForOutput = "C:\\Users\\bernhard.konig\\Documents\\ASync\\home\\Privat\\SAV\\Fachgruppe Data Science\\Paper July 2018\\"
folderForOutput = "H:\\Privat\\SAV\\Fachgruppe Data Science\\Paper July 2018\\"
@assert isdir(folderForOutput) "Directory does not exist: $(folderForOutput)"
# define functions
# include(joinpath(@__DIR__,"..","tutorials","5.ReservingExample_functions.jl"))
# include(joinpath(Pkg.dir("DecisionTrees"),"tutorials","5.ReservingExample_functions.jl")) #this may not use the most current file but rather the file in the folder with compiled julia code, e.g. C:\\Users\\your_username\\.julia\\packages\\DecisionTrees\\3AXAR\\.... or similar
fn_file = joinpath(pwd(), "tutorials", "5.ReservingExample_functions.jl")
@assert isfile(fn_file)
include(fn_file)

##############################
# Read the data
##############################
# note: the underyling data was generated with R (consider the readme.txt in the folder data\\reservingAccident\\)
datafile = joinpath("data", "reservingAccident", "Simulated.Cashflows.csv")
@assert isfile(datafile) "File not found: $(datafile). Consider the readme.txt in the folder data\\reservingAccident\\"
elt = [String,String,String,Int,Int,Int,String,Int]
elt = vcat(elt, repeat([Float64], 12))
@info("Reading data...")
@time fullData = CSV.read(datafile, DataFrame, types=elt, strict=true);

#=
ecdf11=ecdf(fullData[:PayCum11])
using ORCA
using PlotlyJS
x_values=collect(0:1000:maximum(ecdf11))
x_values=collect(0:1000:50000)
y_values=ecdf11.(x_values)
layout = Layout(;title="Cumulative claim count distribution",ylabel="Percentage of Claims",font_size=23)
pl=PlotlyJS.plot(x_values,y_values,layout,xaxis=("Claim Size"),ticks = :native,dpi=200,show=true)
#PlotlyJS.savefig(pl,"c:\\temp\\plot.jpeg")

#consider claim count in percentage versus total aggregate loss amount (in percentage)
#plot cumulative distribution of PayCum11
vsorted=sort(deepcopy(fullData[:PayCum11]))
tbl=count_vs_agg_amount(vsorted,1000)
writedlm(joinpath(folderForOutput,"ExcelTables","agg_claim_distribution.csv"),tbl,',')

paycum11sorted=deepcopy(fullData[:PayCum11])
sort!(paycum11sorted)
vsum=sum(paycum11sorted)
@btime aggLossPctSortedInput(.99,$paycum11sorted,vsum) =#

# for development only: define random subset of data
debug = false
if debug
    @warn("Subsetting data:")
    rnd = rand(size(fullData, 1))
    srtTmp = sortperm(rnd)
    subsetIdx = srtTmp[1:200000]
    fullData = fullData[subsetIdx,:]
    folderForOutput = "C:\\temp\\"
end

# Define a random split into training and validation data
# training shall be 70% of the data
Random.seed!(1239)
trnValCol = BitArray(rand() < 0.7 for x in 1:size(fullData, 1));
fullData[:trnValCol] = trnValCol

# consider deciles (and 'twenty'-iles) of ultimate claim amounts
decilesThresholds = quantile(fullData[:PayCum11], collect(range(0.1, step=.1, stop=.91)))
twentythresholds = quantile(fullData[:PayCum11], collect(range(0.05, step=.05, stop=.96)))

dec = map(x->searchsortedfirst(decilesThresholds, x), fullData[:PayCum11])
twen = map(x->searchsortedfirst(twentythresholds, x), fullData[:PayCum11])
fullData[:tenGroupsSortedByUltimate] = dec
fullData[:twentyGroupsSortedByUltimate] = twen

@warn("We may want to redo the whole calculation with no validation data -> the CL factors will be based on all of the data (not only the training porportion).")

# The data should match this total for the column Paycum11
expectedSum = 9061763108
if !debug
    @assert sum(fullData[:PayCum11]) == expectedSum "Expected $(expectedSum), have: $(sum(fullData[:PayCum11]))"
end
# The data is described in the paper mentioned above

# We provide a short summary from the Readme of Wuethrich and Gabrielli
# The variable ClNr is the claim number, this is a unique numerical identifier for each individual claim.
# LoB is the line of business which should be of factor type (categorical) and may take the labels 1; : : : ; 4.
# cc is the claims code that should be of factor type (categorical) and may take the labels 1; : : : ; 53 (with possible gaps).
# AY is the year of claims occurrence (accident year) that should be of integer type and may take the values 1994; : : : ; 2005.
# AQ is the quarter of claims occurrence (accident quarter) that should be of integer type and may take the values 1; : : : ; 4.
# age is the age of the injured that should be of integer type and may take the values 15; : : : ; 70.
# inj_part is the body part injured that should be of factor type (categorical) and may take the labels 1; : : : ; 99 (with gaps).

# see what Plot backend is used
# backend() #e.g. Plots.GRBackend

# set independent variables
selected_explanatory_vars = ["AQ","LoB","age","cc","inj_part"]
categoricalVars = ["LoB","inj_part","cc"]

# check type of each column
for x in vcat(selected_explanatory_vars, "AY")
    println(x)
    println(eltype(fullData[Symbol(x)]))
    println("first ten values")
    print(fullData[1:10,Symbol(x)])

    # count frequencies
    ft = countmap(fullData[Symbol(x)])
    # sort data to be plotted
    ks = collect(keys(ft))
    srt = sortperm(ks)
    vs = collect(values(ft))
    ks = ks[srt]
    vs = vs[srt]
    # plot the frequency count
    StatsPlots.bar(ks, vs, title="", label="", xaxis=x, yaxis="number of claims", fmt=:svg)
    savedir = joinpath(folderForOutput, "figures") # folderForOutput #"C:\\temp\\"
    isdir(savedir) && StatsPlots.savefig(joinpath(savedir, string(x, "_graph.pdf")))
    println("")
end

# consider unique elements for each variable
for x in 1:size(fullData, 2)
    println(x)
    println(names(fullData)[x])
    @show size(unique(fullData[:,x]))
end

#=
#consider a sample of the data (PayCum01-PayCum11)
fullData[1:400,10:end-1]

any(fullData[:PayCum01].!=fullData[:PayCum02])
any(fullData[:PayCum10].!=fullData[:PayCum11]) =#

#=
    #Let us see whether all claims are 'increasing' in time (i.e. whether all incremental payments are non-negative)
    #NOTE: This may take a few minutes to compute!
    #->there are some rows/claims which have non increasing payments over time
    #99.7% of the claims are nondecreasing though

    isNonDecreasing=checkIfPaymentsAreIncreasing(fullData[1:30])
    sum(isNonDecreasing)/length(isNonDecreasing) # 99.7% of the claims are nondecreasing
    nonDecreasingRows=.!isNonDecreasing
    nonDecreasingRowsInteger=find(nonDecreasingRows)
    @show size(nonDecreasingRowsInteger)
    #claims which are decreasing over time (at some point):
    #fullData[nonDecreasingRowsInteger,:] =#

##############################
# Prepare the data
##############################
# consider paid to date amount
paidToDateFullDataPerRow = getPaidToDatePerRow(fullData)
fullData[:paidToDate] = paidToDateFullDataPerRow

# consider the Accident Years
ayears = sort(unique(fullData[:AY]))
maxAY = maximum(ayears)
@assert maxAY == 2005

# prepare the whole dataset in the same manner as the training data set
    # indicator index whether claim is known by year end 2005
indexOfClaimIsReportedByYE2005 = fullData[:AY] + fullData[:RepDel] .<= 2005
dataKnownByYE2005 = fullData[indexOfClaimIsReportedByYE2005,:]
    # paycumcolumns=[Symbol(string("PayCum",lpad(ii,2,0))) for ii=0:11]
    # map(x->sum(fullData[.!indexOfClaimIsReportedByYE2005,x]),paycumcolumns)
    # aggregate(df,:someColumn,[sum, mean])
dtmKnownByYE2005, sett, dfpreppedUnused = prepare_dataframe_for_dtm!(dataKnownByYE2005, keycol="ClNr", numcol="PayCum00", trnvalcol="trnValCol", independent_vars=selected_explanatory_vars, treat_as_categorical_variable=categoricalVars);
    # define general model settings
updateSettingsMod!(sett, ignoreZeroDenominatorValues=true, model_type="build_tree", write_dot_graph=true, writeTree=false, graphvizexecutable="C:\\Program Files (x86)\\Graphviz2.38\\bin\\dot.exe")
    # a note on the data
    # dataKnownByYE2005 is a DataFrame (it contains more data such as all PayCum columns and the AY)
    # dtmKnownByYE2005 is a DTMTable struct (it has a different structre than the DataFrame)

# consider paidToDate and 'truth'
paidToDatePerRow = getPaidToDatePerRow(dataKnownByYE2005)
# trueUltimatePerRow (this is not meaningful because we have two sets of claims: fullData[:PayCum11] AND dtmKnownByYE2005[:PayCum11])
# trueReservePerRow=trueUltimatePerRow.-paidToDatePerRow #this is not meaningful either

# build 'triangle' (in this case, we will actually get a rectangle rather than a triangle)
@info("Fitting 'plain' CL to all data (including 'unknown' part of the triangle)")
triangleInclIBNyRClaims = buildTriangle(fullData)
trueUltimate = triangleInclIBNyRClaims[:,end]

# CL
# consider aggregate CL Model
@info("Fitting 'plain' CL to all data")
triangle = buildTriangle(dataKnownByYE2005)
clAllLOBs = chainLadder(triangle)

trueReserves = trueUltimate .- clAllLOBs[:paidToDate]

# Consider CL Ultimate per known claim
# LDF = Loss Development Factor (also chain ladder factor)
ayperRow = dataKnownByYE2005[:AY]
CLfactorsToUltimate = clAllLOBs[:factorsToUltimate]
CLLDFperRow = map(x->CLfactorsToUltimate[x - 1993], ayperRow)
CLUltimatePerRow = CLLDFperRow .* paidToDatePerRow
cldf = DataFrame(ldf=CLLDFperRow, ultimate=CLUltimatePerRow, paidToDate=paidToDatePerRow)
cldf[:reserves] = cldf[:ultimate] .- cldf[:paidToDate]
CLcomparisons = CL_est_per_variable(fullData, cldf)
summaryByVar = vcat(collect(values(CLcomparisons)))
clResultsFolder = joinpath(folderForOutput, "clResults")
isdir(clResultsFolder) && CSV.write(joinpath(clResultsFolder, "resultsCL.csv"), summaryByVar)

if !debug
    @assert isapprox(0, sum(CLUltimatePerRow) - sum(clAllLOBs[:ultimate]), atol=1e-4) # should be zero
    for i in ayears
        idx = ayperRow .== i
        @assert isapprox(0, sum(CLUltimatePerRow[idx]) - clAllLOBs[:ultimate][i - 1993], atol=1e-4) # should be zero for each year
    end
end
CLTotalUltimate = sum(CLUltimatePerRow)

# consider CL per LoB
@info("Fitting 'plain' CL models to each LoB")
LOB = "1";triangleLOB1 = buildTriangle(dataKnownByYE2005[dataKnownByYE2005[:LoB] .== LOB,:]);clLOB1 = chainLadder(triangleLOB1);trueUltimateLOB1 = buildTriangle(fullData[fullData[:LoB] .== LOB,:])[:,end];clLOB1[:trueReserves] = trueUltimateLOB1 .- clLOB1[:paidToDate];clLOB1[:error] = clLOB1[:trueReserves] .- clLOB1[:reserves];
LOB = "2";triangleLOB2 = buildTriangle(dataKnownByYE2005[dataKnownByYE2005[:LoB] .== LOB,:]);clLOB2 = chainLadder(triangleLOB2);trueUltimateLOB2 = buildTriangle(fullData[fullData[:LoB] .== LOB,:])[:,end];clLOB2[:trueReserves] = trueUltimateLOB2 .- clLOB2[:paidToDate];clLOB2[:error] = clLOB2[:trueReserves] .- clLOB2[:reserves];
LOB = "3";triangleLOB3 = buildTriangle(dataKnownByYE2005[dataKnownByYE2005[:LoB] .== LOB,:]);clLOB3 = chainLadder(triangleLOB3);trueUltimateLOB3 = buildTriangle(fullData[fullData[:LoB] .== LOB,:])[:,end];clLOB3[:trueReserves] = trueUltimateLOB3 .- clLOB3[:paidToDate];clLOB3[:error] = clLOB3[:trueReserves] .- clLOB3[:reserves];
LOB = "4";triangleLOB4 = buildTriangle(dataKnownByYE2005[dataKnownByYE2005[:LoB] .== LOB,:]);clLOB4 = chainLadder(triangleLOB4);trueUltimateLOB4 = buildTriangle(fullData[fullData[:LoB] .== LOB,:])[:,end];clLOB4[:trueReserves] = trueUltimateLOB4 .- clLOB4[:paidToDate];clLOB4[:error] = clLOB4[:trueReserves] .- clLOB4[:reserves];

addTotal!(clLOB1)
addTotal!(clLOB2)
addTotal!(clLOB3)
addTotal!(clLOB4)

trueReservesPerLOB = Dict{String,Array{Float64}}()
trueReservesPerLOB["1"] = clLOB1[:trueReserves]
trueReservesPerLOB["2"] = clLOB2[:trueReserves]
trueReservesPerLOB["3"] = clLOB3[:trueReserves]
trueReservesPerLOB["4"] = clLOB4[:trueReserves]

clPerLOB = Dict{String,DataFrame}()
clPerLOB["1"] = clLOB1
clPerLOB["2"] = clLOB2
clPerLOB["3"] = clLOB3
clPerLOB["4"] = clLOB4

#=
for i=1:4;tmp=@eval $(Symbol(string("triangleLOB",i)));writedlm(joinpath(folderForOutput,"Arius","triangle_LOB$(i).csv"),tmp,',');end; =#

##
# Combine the four CL models on the whole data
# LDF = Loss Development Factor (also chain ladder factor)
ayAndLOBperRow = dataKnownByYE2005[[:AY,:LoB]]
CLLDFperRow_byLOB = map(x->clPerLOB[String(ayAndLOBperRow[x,2])][:factorsToUltimate][ayAndLOBperRow[x,1] - 1993], 1:size(ayAndLOBperRow, 1))
CLUltimatePerRow_byLOB = CLLDFperRow_byLOB .* paidToDatePerRow
cldf_byLOB = DataFrame(ldf=CLLDFperRow_byLOB, ultimate=CLUltimatePerRow_byLOB, paidToDate=paidToDatePerRow)
cldf_byLOB[:reserves] = cldf_byLOB[:ultimate] .- cldf_byLOB[:paidToDate]
CLcomparisons_byLOB = CL_est_per_variable(fullData, cldf_byLOB)
summaryByVar_byLOB = vcat(collect(values(CLcomparisons_byLOB)))
isdir(clResultsFolder) && CSV.write(joinpath(clResultsFolder, "resultsCL_by_LOB.csv"), summaryByVar_byLOB)

##

if !debug
    @assert all(isapprox.(vcat(trueReserves, sum(trueReserves)) .- clLOB1[:trueReserves] .- clLOB2[:trueReserves] .- clLOB3[:trueReserves] .- clLOB4[:trueReserves], 0))
    @assert all(isapprox.(trueUltimate .- trueUltimateLOB1 .- trueUltimateLOB2 .- trueUltimateLOB3 .- trueUltimateLOB4, 0))
end

##############################
# Consider CL decision trees
##############################

# Define minimum weight per leaf (depending on the development year)
# note:these weights are calibrated for a total training weight of 3477583 (Notably the number of claims will be smaller for higher LDFs!)
@show sum(dataKnownByYE2005[:trnValCol]) # expected value: 3477583
@assert debug || (sum(dataKnownByYE2005[:trnValCol]) == 3477583)

# each entry corresponds to an LDF Year
modelsWeightsPerLDF = Vector{Vector{Float64}}()
# the values were selected by expert judgement
# primarily we considered reversals (i.e when the observed ratio in the validation data is not monotonic in the number of leaves (as the leaf number is defined such that the training observed ratios are increasing))
push!(modelsWeightsPerLDF,[150000,200000,270000,250000,235000,210000,200000,240000,180000,130000,130000])
# modelsWeightsPerLDF=Vector{Vector{Float64}}()
# alternative
# push!(modelsWeightsPerLDF,[150000,200000,270000,220000,235000,210000,200000,200000,130000,130000,15000000])

# alternatives

# old weights for sse
# push!(modelsWeightsPerLDF,[130000,170000,240000,220000,235000,210000,200000,240000,130000,130000,15000000])

# You can provide additional sets of minimum weights to see the differences between two models
# push!(modelsWeightsPerLDF,[70000,170000,200000,220000,235000,210000,200000,240000,200000,130000,15000000])

# WEIGHTS FOR DIFFERENCE FN: push!(modelsWeightsPerLDF,[130000,170000,240000,220000,235000,210000,200000,240000,130000,130000,15000000])
# version where the weight is increasing
# push!(modelsWeightsPerLDF,[80000,170000,240000,220000,240000,240000,240000,240000,240000,240000,15000000])

for i = 1:length(modelsWeightsPerLDF)
    @assert length(modelsWeightsPerLDF[i]) == length(ayears) - 1 # == number of LDFs
end

# goal: build a tree for each LDF (i.e. a tree for 12-24 months, a second tree for 24-36 months, etc.)

treeResults = Dict{Float64,DataFrame}()
treeResultsAgg = Dict{Float64,DataFrame}()

@info("Starting model fit...")
@show size(modelsWeightsPerLDF)
# folderForOutput should be defined further above
folderForSingleTreeResults = joinpath(folderForOutput, "figures", "SingleTreeApproach")
@assert isdir(folderForSingleTreeResults) "Directory does not exist: $(folderForSingleTreeResults)"
# actual model run (this may take a few minutes)
updateSettingsMod!(sett,crit="sse",model_type="build_tree",subsampling_features_prop=1.0)
@time this_ldfarr = runModels!(dataKnownByYE2005, dtmKnownByYE2005, modelsWeightsPerLDF, treeResults, treeResultsAgg, selected_explanatory_vars, categoricalVars, folderForSingleTreeResults, clAllLOBs, paidToDatePerRow, sett);

@info("Finished building models.")
# write data to disk
mdl_stats = write_results(treeResults, treeResultsAgg, folderForSingleTreeResults, trueReservesPerLOB, clPerLOB)

############################################################
# boosted model
############################################################
@info("Starting Boosting models...")
folderForBoostingResults = joinpath(folderForOutput, "figures", "boostingResults")
@assert isdir(folderForBoostingResults) "Directory does not exist: $(folderForBoostingResults)"
treeResultsBoosting = Dict{Float64,DataFrame}()
treeResultsAggBoosting = Dict{Float64,DataFrame}()
updateSettingsMod!(sett,crit="sse",iterations=8,learningRate=.15,model_type="boosted_tree",subsampling_features_prop=0.6)
@time ldfArrBoosting = runModels!(dataKnownByYE2005, dtmKnownByYE2005, modelsWeightsPerLDF, treeResultsBoosting, treeResultsAggBoosting, selected_explanatory_vars, categoricalVars, folderForBoostingResults, clAllLOBs, paidToDatePerRow, sett);

mdl_statsBoosting = write_results(treeResultsBoosting, treeResultsAggBoosting, folderForBoostingResults, trueReservesPerLOB, clPerLOB)
# Note that the boosing model results in three different estimates
# the dictionary mdl_statsBoosting has entries 1.0, 1.2 and 1.4.
# 1.0 canonical estimate from the boosting model
# 1.4 the estimates were grouped in N groups of equal size (N is user parameter, default is 1000) and then averaged, thus there are at most 1000 distinct estimated values
# 1.2 a moving average is applied to the fitted values of "1.4" above

@warn("If we were to update the 'tree-CL' factors such that they are based on all data (currently they are only using the training data), the model might further improve.")

#=
#try some alternative models
NN=20
this_ldfarr[1:NN,1]
ldfArr[1:NN,1]
    ldfYear=1
    selectedWeight=150000
    settC=deepcopy(sett)
    updateSettingsMod!(settC,crit="sse",subsampling_features_prop=1.0)
    resultingFiles,resM= runSingleModel(dataKnownByYE2005,dtmKnownByYE2005,selectedWeight,ldfYear,2005,               selected_explanatory_vars,categoricalVars,"C:\\temp\\331\\",settC)
    0

#tbd: compare mse versus sse
selectedWeight=270000
while (selectedWeight < 180000)
    ldfYear=1
    #selectedWeight=150000
    settC=deepcopy(sett)
    updateSettingsMod!(settC,crit="mse")
    resultingFiles,resM= runSingleModel(dataKnownByYE2005,dtmKnownByYE2005,selectedWeight,ldfYear,2005,               selected_explanatory_vars,categoricalVars,"C:\\temp\\331\\",settC)

    settC=deepcopy(sett)
    updateSettingsMod!(settC,crit="sse",subsampling_features_prop=1.0)
    resultingFiles,resM= runSingleModel(dataKnownByYE2005,dtmKnownByYE2005,selectedWeight,ldfYear,2005,               selected_explanatory_vars,categoricalVars,"C:\\temp\\331\\",settC)
    selectedWeight+=10000
end

ldfYear=9
selectedWeight=180000
settC=deepcopy(sett)
updateSettingsMod!(settC,crit="difference",iterations=8,learningRate=.15,model_type="boosted_tree",subsampling_features_prop=0.6)
resultingFiles,resM= runSingleModel(dataKnownByYE2005,dtmKnownByYE2005,selectedWeight,ldfYear,2005,               selected_explanatory_vars,categoricalVars,"C:\\temp\\331\\",settC)
dfpred=predict(resM,dtmKnownByYE2005.features)
fittedValues=dfpred[:RawEstimate]
@show mean(fittedValues),median(fittedValues)
0

ldfYear=1
selectedWeight=130000
settC=deepcopy(sett)
updateSettingsMod!(settC,crit="gamma")
resultingFiles,resM= runSingleModel(dataKnownByYE2005,dtmKnownByYE2005,selectedWeight,ldfYear,2005,               selected_explanatory_vars,categoricalVars,"C:\\temp\\331\\",settC)
0 =#

#=
#test algorithm on smaller data Set

if debug
    @warn("Subsetting data:")
    this_ldfy=1
    this_selectedWeight=Int(round(1/size(fullData,1)*modelsWeightsPerLDF[1][this_ldfy]*trnsize,digits=0))
    this_selectedWeight=5500
    trnsize=40000
    cvsampler=CVOptions(0,0.0,true)
    test_on_subset(dataKnownByYE2005,this_ldfy,folderForOutput,sett,trnsize,this_selectedWeight,cvsampler)

    cvsampler=CVOptions(-2,0.0,true)
    test_on_subset(dataKnownByYE2005,this_ldfy,folderForOutput,sett,trnsize,this_selectedWeight,cvsampler)
0
end =#

#
# Alternative Model: build a (plain) CL model for each Leaf indicated by the decision tree for the first develoment year
#

folderForOneCLModelPerLeaf = joinpath(splitdir(folderForSingleTreeResults)[1], "oneCLModelPerLeaf")
@assert isdir(folderForOneCLModelPerLeaf)

ldfYear = 1
selectedWeight = 150000
settC = deepcopy(sett)
updateSettingsMod!(settC,crit="sse",model_type="build_tree",subsampling_features_prop=1.0)
resultingFiles, resM = runSingleModel(dataKnownByYE2005, dtmKnownByYE2005, selectedWeight, ldfYear, 2005, selected_explanatory_vars, categoricalVars, folderForOneCLModelPerLeaf, settC)

fitted, leafNrs = DecisionTrees.predict(resM, dtmKnownByYE2005.features)
list_of_leafnumbers = sort(unique(leafNrs))
thisLeafNr = list_of_leafnumbers[3]
clPerLeaf = Dict{Int,DataFrame}()
for thisLeafNr in list_of_leafnumbers
    idx = leafNrs .== thisLeafNr
    subsetOfData = dataKnownByYE2005[idx,:]
    subTriangle = buildTriangle(subsetOfData)
    # write triangle to disk
    writedlm(joinpath(folderForOneCLModelPerLeaf, "triangle_of_leaf_number_$(thisLeafNr).csv"), subTriangle, ',')
    subCL = chainLadder(subTriangle)
    subTrueUltimate = subTriangle[:,end]
    subCL[:trueReserves] = subTrueUltimate .- subCL[:paidToDate]
    subCL[:error] = subCL[:trueReserves] .- subCL[:reserves]
    # fill dictionary
    clPerLeaf[thisLeafNr] = deepcopy(subCL)
end

# derive statistics
ayAndLeafNrPerRow = hcat(dataKnownByYE2005[:AY], leafNrs)
CLLDFperRow_byLeaf = map(x->clPerLeaf[ayAndLeafNrPerRow[x,2]][:factorsToUltimate][ayAndLeafNrPerRow[x,1] - 1993], 1:size(ayAndLeafNrPerRow, 1))
CLUltimatePerRow_byLeaf = CLLDFperRow_byLeaf .* paidToDatePerRow
cldf_byLeaf = DataFrame(ldf=CLLDFperRow_byLeaf, ultimate=CLUltimatePerRow_byLeaf, paidToDate=paidToDatePerRow)
cldf_byLeaf[:reserves] = cldf_byLeaf[:ultimate] .- cldf_byLeaf[:paidToDate]
CLcomparisons_byLeaf = CL_est_per_variable(fullData, cldf_byLeaf)
summaryByVar_byLeaf = vcat(collect(values(CLcomparisons_byLeaf)))
isdir(folderForOneCLModelPerLeaf) && CSV.write(joinpath(folderForOneCLModelPerLeaf, "resultsCL_by_LeafNr.csv"), summaryByVar_byLeaf)


text:
# Tests for mutual information function
# (NB These tests will fail with the shrinkage estimator, because the 1D,
# 2D and 3D distributions on the right will shrink towards different
# targets by default, whereas the mutual information formula overcomes
# by estimating the probabilities of the 3D distribution, then projecting)

# Test mutual information formula
@test get_mutual_information(arr1, arr2) ≈ get_entropy(arr1) + get_entropy(arr2) - get_entropy(arr1, arr2)
println("Mutual information formula passed.")


text:
InputFolder = './Images/Yeasts/';
OutputFolder = './Results/Images/Yeasts3/';
Rescale = 2;

@iA = '*.tif';

@fxg_mMultiPassThr [iA] > [M];
params.NeighRad = 25;
params.GaussRad = 1;
params.Rescale = 1;
params.Th1 = 6;
params.Th2 = 10;
params.Rel = 0.9;
params.MinArea = 200;
/endf

@fxm_mModBinWat [M] > [M];
params.SmallHolesArea = 0;
params.DistStdRad = 7;
params.MinDistLocVar = 9;
/endf

/show iA > M;
/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
export TRMA97_abs

function TRMA97_abs(nlp 		:: AbstractNLPModel,
                 	nlpstop 	:: NLPStopping;
				    kwargs...
               		)

	return TRARC2(nlp, nlpstop; TR = TrustRegion(10.0),
				  c = Combi(hessian_dense, PDataMA97{eltype(nlp.meta.x0)}, solve_modelTRDiagAbs, preprocessMA97, decreaseFact, Tparam{eltype(nlp.meta.x0)}()),
				  kwargs...
				  )
end


text:
# This script produces all of the plots that appear in the README.md file
#
# Requires Plots package
using Random
using Distributions, Plots
using Changepoints

Random.seed!(1) # Fix random seed

#########################
# Simulate changepoints #
#########################

# Simulate time series with changepoints
n = 1000          # Sample size
λ = 70            # Freq of changepoints

# Mu parameter sampled from a N(0, 1) distribution for each changepoint
μ, σ = Normal(0,1), 1.0
data, cps = @changepoint_sampler n λ Normal(μ, σ)
changepoint_plot(data, cps)
savefig("Plots_example.png")

# PELT

pelt_cps, cost = @PELT data Normal(:?, 1.0)
changepoint_plot(data, pelt_cps)
savefig("Plots_example_pelt.png")

# BS

bs_cps = @BS data Normal(:?, 1.0)
changepoint_plot(data, bs_cps[1])
savefig("Plots_example_BS.png")

# CROPS

seg_cost = NormalMeanSegment(data, σ)
pen1, pen2 = 4.0, 100.0 # Penalty range
crops_output = CROPS(seg_cost, n, pen1, pen2)
elbow_plot(crops_output)
savefig("Plots_elbow_plot.png")

# MOSUM

G = 35
MOSUM_output = MOSUM(data, G)
mosum_plot(MOSUM_output)
savefig("Plots_mosum_plot.png")

# MOSUM multi-scale

Gset = [20, 30, 50, 80, 130]
MOSUM_multi_scale_output =  MOSUM_multi_scale(data, Gset)
changepoint_plot(data, MOSUM_multi_scale_output)
savefig("Plots_mosum_multi_scale.png")

# WBS
WBS_return = WBS(data, min_length = 10 )
#WBS_return = Changepoints.@WBS data
seg_cost_sSIC = Changepoints.sSIC(data)
WBS_cps = get_WBS_changepoints(seg_cost_sSIC, WBS_return, 14)
changepoint_plot(data, WBS_cps[1])
savefig("Plots_WBS.png")

SeedBS_return = WBS(data, do_seeded = true, min_length = 10)
#SeedBS_return = @WBS data do_seeded=true
SeedBS_cps = get_WBS_changepoints(seg_cost_sSIC, SeedBS_return, 14)
changepoint_plot(data, SeedBS_cps[1])
savefig("Plots_SeedBS.png")


text:
construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{<:PSY.ACBranch, <:AbstractBranchFormulation},
    ::Union{Type{CopperPlatePowerModel}, Type{AreaBalancePowerModel}},
) = nothing

construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{<:PSY.DCBranch, <:AbstractDCLineFormulation},
    ::Union{Type{CopperPlatePowerModel}, Type{AreaBalancePowerModel}},
) = nothing

# This method might be redundant but added for completness of the formulations
construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::DeviceModel{<:PSY.Branch, <:Type{UnboundedBranches}},
    ::Type{<:PM.AbstractPowerModel},
) = nothing

# For DC Power only. Implements Bounds only and constraints
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{B, Br},
    ::Type{S},
) where {
    B <: PSY.ACBranch,
    Br <: AbstractBoundedBranchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end
    !isnothing(get_feedforward(model)) &&
        throw(IS.ConflictingInputsError("$(Br) formulation doesn't support FeedForward. Use Constrained Branch Formulation instead"))
    branch_rate_bounds!(psi_container, devices, model, S)
    branch_rate_constraints!(psi_container, devices, model, S, get_feedforward(model))
    return
end

# For DC Power only. Implements Constraints only
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end
    branch_rate_constraints!(psi_container, devices, model, S, get_feedforward(model))
    return
end

# For AC Power only. Implements Bounds on the active power and rating constraints on the aparent power
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{B, <:AbstractBranchFormulation},
    ::Type{S},
) where {B <: PSY.ACBranch, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end
    branch_rate_bounds!(psi_container, devices, model, S)
    branch_rate_constraints!(psi_container, devices, model, S, get_feedforward(model))
    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{B, Br},
    ::Type{S},
) where {B <: PSY.DCBranch, Br <: AbstractDCLineFormulation, S <: PM.AbstractPowerModel}
    devices = get_available_components(B, sys)
    if !validate_available_devices(B, devices)
        return
    end
    branch_rate_constraints!(psi_container, devices, model, S, get_feedforward(model))
    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{PSY.MonitoredLine, FlowMonitoredLine},
    ::Type{S},
) where {S <: PM.AbstractActivePowerModel}
    devices = get_available_components(PSY.MonitoredLine, sys)
    if !validate_available_devices(PSY.MonitoredLine, devices)
        return
    end
    branch_flow_constraints!(psi_container, devices, model, S, get_feedforward(model))
    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{PSY.MonitoredLine, FlowMonitoredLine},
    ::Type{S},
) where {S <: PM.AbstractPowerModel}
    devices = get_available_components(PSY.MonitoredLine, sys)
    if !validate_available_devices(PSY.MonitoredLine, devices)
        return
    end
    branch_rate_constraints!(psi_container, devices, model, S, get_feedforward(model))
    branch_flow_constraints!(psi_container, devices, model, S, get_feedforward(model))
    return
end


text:
"""
Data module for preparing data for analysis with
    - GaussianProcess (src/gpr/GaussianProcess.jl)
"""

module Data

using OceanTurb
using OrderedCollections

include("../LES/custom_avg.jl")
export custom_avg

# harvesting Oceananigans data
include("../LES/read_les_output.jl")
export get_les_data

# running OceanTurb KPP simulations based on OceananigansData conditions
include("../KPP/run.jl")
export closure_free_convection_kpp_full_evolution,
       closure_free_convection_kpp

# running OceanTurb TKE simulations based on OceananigansData conditions
include("../TKE/run.jl")
export closure_free_convection_tke_full_evolution,
       closure_free_convection_tke

# normalization
include("scalings.jl")
export  min_max_scaling,
        scale, # normalize
        unscale # un-normalize

# pre- and post-processing on the normalized data
include("problems.jl")
export  Problem,
        Sequential,
        Residual,
        Slack,
        SequentialProblem,
        ResidualProblem,
        SlackProblem,
        get_problem_v,
        predict # important!

export Sequential_KPP,
       Sequential_TKE

include("sequential.jl")
include("residual.jl")
include("slack.jl")
export  get_predictors_targets,
        postprocess_prediction

include("modify_predictor_fns.jl")
export append_tke,
       partial_temp_profile

include("convective_adjust.jl")
export convective_adjust!

# ProfileData struct
export  ProfileData,
        data

function approx_initial_buoyancy_stratification(T_initial,z)
    # b_initial = 𝒟.T[:,1] .* α*g
    α = 2e-4
    g = 9.80665
    N² = (T_initial[1] - T_initial[end])*α*g / z[1] # approximate initial buoyancy gradient N², where b = N²z + 20*α*g and T = N²z/(αg) + 20
    return N²
end

struct StateVariables
    tke_avg # turbulent kinetic energy
end

"""
ProfileData
# Description
- data structure for preparing profile data for analysis with gpr or nn.
# Data Structure and Description
    v::Array,           Nz x Nt array of T or wT values directly from the LES simulation, not preprocessed.
    vavg::Array,        Nt-length array of Nz-length vectors from LES simulation, scaled and pre-processed
    x::Array,           all simulation inputs, scaled and pre-processed
    y::Array,           all simulation inputs, scaled and pre-processed
    x_train::Array,     training inputs (predictors; array of states). (length-n array of D-length vectors, where D is the length of each input n is the number of training points)
    y_train::Array,     training outputs (predictions) (length-n array of D-length vectors).
    validation_set::Array, vector of indices corresponding to validation data
    z::Array,           Nz-length vector of depth values
    zavg::Array,        length-D vector; depth values averaged to D gridpoints
    t::Array,           timeseries [seconds]
    Nt::Int64,          length(timeseries)
    problem::Problem,   what mapping you wish to evaluate with the model. (Sequential("T"), Sequential("wT"), Residual("T"), Residual("KPP"), or Residual("TKE"))
    all_problems::Array{Array{Any,1},1}
    state_variables::StateVariables
    modify_predictor_fn::Function
    convective_adjust::Function
"""
struct ProfileData
    v       ::Array{Float64,2}
    vavg    ::Array{Array{Float64,1},1}
    x       ::Array{Array{Float64,1},1}
    y       ::Array{Array{Float64,1},1}
    x_train ::Array{Array{Float64,1},1}
    y_train ::Array{Array{Float64,1},1}
    validation_set::Vector{Int64}
    z       ::Vector{Float64}
    zavg    ::Vector{Float64}
    t       ::Vector{Float64}
    Nt      ::Int64
    problem ::Problem
    all_problems::Array{Array{Any,1},1}
    state_variables::StateVariables
    modify_predictor_fn::Function
    convective_adjust::Function
end

"""
data(filename, problem; D=16, N=4)

# Description
Returns a ProfileData object based on data from `filename`

# Arguments
- `filename`: (string)  Name of the NetCDF or JLD2 file containing the data from the Oceananigans simulation.
- `problem`: (Problem). What mapping you wish to evaluate with the model.

# Keyword Arguments
- `D`: (integer)        Number of gridpoints in the z direction to average the data to for training and prediction.
- `N`: (integer)        Interval between the timesteps to be reserved for training data (default 4).
                        If N=4, the profile data for every 4 timesteps will be reserved for training (~25% training data);
                        the rest will be used in the validation set.
"""
function data(filename::String, problem::Problem; D=16, N=4)

    # collect data from Oceananigans simulation output file
    les = get_les_data(filename) # <: OceananigansData

    # eddy diffusivity
    κₑ = les.κₑ

    # depth values
    z = les.z
    zavg = custom_avg(z, D)

    # get v (variable array, Nz x Nt) and cut out the first 2 hours
    start = floor(Int64, 7200 / (les.t[2] - les.t[1]))+1
    s = start:length(les.t)
    v = get_v(problem, les)[:,s]

    # timeseries [s]
    t = les.t[s]
    Nt = length(t)

    # approximate buoyancy stratification at the initial timestep
    N² = approx_initial_buoyancy_stratification(les.T[:,1],z)

    # state variables
    array_to_avg(array) = [custom_avg(array[:,j], D) for j in s]
    state_variables = StateVariables(array_to_avg(les.tke))

    # modify_predictor_fn
    modify_predictor_fn(state, time_index) = problem.modify_predictor_fn(state, time_index, state_variables)

    # convective_adjust
    if problem.convective_adjust
        convective_adjust = convective_adjust!
    else
        convective_adjust = x -> x
    end

    # get problem (sets how the data will be pre- and post-processed)
    specific_problem = get_problem(problem, les, v, N², D, Nt)

    # compress variable array to D gridpoints to get an Nt-length array of D-length vectors
    vavg = [custom_avg(v[:,j], D) for j in 1:Nt]

    # divide up the data
    total_set = 1:(Nt-1)
    training_set = 1:N:(Nt-1)
    validation_set = setdiff(total_set, training_set)
    n_train = length(training_set)

    # get the predictors and target predictions to train the model on
    x, y = get_predictors_targets(vavg, specific_problem)

    x_train = [modify_predictor_fn(x[i], i) for i in training_set]
    y_train = y[training_set]

    # if the data contains multiple files, we'll need a way to postprocess the data separately for each one
    # (preprocessing is already taken care of before the data is merged)
    all_problems = [[specific_problem, length(x)]]

    return ProfileData(v, vavg, x, y, x_train, y_train, validation_set, z, zavg, t, Nt, specific_problem, all_problems, state_variables, modify_predictor_fn, convective_adjust)
end

"""
data(filename, D; N=4)

# Description
Returns an instance of ProfileData containing training data from multiple simulations.
*** Important:
    ONLY v, x_train and y_train contain data from all filenames;
    the remaining attributes are taken from the first filename in filenames

# Arguments
- `filenames`: (string)  Vector of filenames (.nc or .jld2) to collect data from.
- `problem`: (Problem). What mapping you wish to evaluate with the model. (Sequential("T"), Sequential("wT"), Residual("T"), Residual("KPP"), or Residual("TKE"))

# Keyword Arguments
- `D`: (integer)        Number of gridpoints in the z direction to average the data to for training and prediction.
- `N`: (integer)        Interval between the timesteps to be reserved for training data (default 4).
                        If N=4, the profile data for every 4 timesteps will be reserved for training (~25% training data);
                        the rest will be used in the validation set.
"""
function data(filenames::Vector{String}, problem::Problem; D=16, N=4)

    # combines data from multiple files
    𝒟 = data(filenames[1], problem; D=D, N=N)

    v = 𝒟.v
    vavg = 𝒟.vavg
    x = 𝒟.x
    y = 𝒟.y
    x_train = 𝒟.x_train
    y_train = 𝒟.y_train
    tke_avg = 𝒟.state_variables.tke_avg
    validation_set = 𝒟.validation_set
    all_problems = 𝒟.all_problems

    t = 𝒟.t
    Nt = 𝒟.Nt

    for filename in filenames[2:end]
        𝒟2 = data(filename, problem; D=D, N=N)

        v = hcat(v, 𝒟2.v) # unscaled
        vavg = vcat(vavg, 𝒟2.vavg)
        x = vcat(x, 𝒟2.x)
        y = vcat(y, 𝒟2.y)
        validation_set = vcat(validation_set, 𝒟2.validation_set)
        x_train = vcat(x_train, 𝒟2.x_train)
        y_train = vcat(y_train, 𝒟2.y_train)
        tke_avg = vcat(tke_avg, 𝒟2.state_variables.tke_avg)
        append!(all_problems, 𝒟2.all_problems)
        t = vcat(t, 𝒟2.t)
        Nt += 𝒟2.Nt
    end

    state_variables = StateVariables(tke_avg)
    modify_predictor_fn(state, time_index) = problem.modify_predictor_fn(state, time_index, state_variables)

    # Note the problem is that from the first file in filenames. This is only included so that the problem type can be determined easily.
    return ProfileData(v, vavg, x, y, x_train, y_train, validation_set, 𝒟.z, 𝒟.zavg, t, Nt, 𝒟.problem, all_problems, state_variables, modify_predictor_fn, 𝒟.convective_adjust)
end

end # module

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
module TestMultiScheduler

import Schedulers
using Schedulers.Internal: @record
using Test

using ..Utils: fetch_finished

function sets2()
    if Threads.nthreads() == 1
        workers1 = workers2 = 1:1
    else
        n1 = Threads.nthreads() ÷ 2
        n2 = Threads.nthreads() - n1
        workers1 = 1:n1
        workers2 = (1:n2) .+ n1
    end
    schedulers = [Schedulers.workstealing, Schedulers.prioritized]
    return Iterators.map(Iterators.product(schedulers, schedulers)) do (f1, f2)
        s1 = f1(workers = workers1)
        s2 = f2(workers = workers2)
        (; label = "$f1-$f2", s1, s2)
    end
end

function check_combination2(s1, s2)
    local t1, t2, t3
    @sync begin
        s1ready = Base.Event()
        t1ready = Base.Event()
        @async begin
            Schedulers.open(s1) do s1
                notify(s1ready)
                @record(:s1ready)
                @record(:wait_t1ready_begin)
                wait(t1ready)
                @record(:wait_t1ready_end)
                @record(:s1_waits_t1_begin)
                wait(t1)  # waiting a s1 task in a s1 task
                @record(:s1_waits_t1_end)
                @record(:s1_waits_t3_begin)
                wait(t3)  # waiting a s1 task in a s1 task
                @record(:s1_waits_t3_end)
            end
        end
        Schedulers.open(s2) do s2
            @record(:wait_s1ready_begin)
            wait(s1ready)
            @record(:wait_s1ready_end)
            t1 = Schedulers.spawn(s1) do
                t2 = Schedulers.spawn(s2) do
                    222
                end
                t3 = Schedulers.spawn(s1) do
                    @record(:t3_waits_t2_begin)
                    local y2 = fetch(t2)  # waiting a s2 task in a s1 task
                    @record(:t3_waits_t2_end)
                    y2 + 111
                end
                @record(:t1_waits_t3_begin)
                local y3 = fetch(t3)  # waiting a s1 task in a s1 task
                @record(:t1_waits_t3_end)
                @record(:t1_waits_t2_begin)
                local y2 = fetch(t2)  # waiting a s2 task in a s1 task
                @record(:t1_waits_t2_end)
                y3 - y2
            end
            notify(t1ready)
            @record(:t1ready)
            @record(:t2_waits_t1_begin)
            wait(t1)  # waiting a s2 task in a s1 task
            @record(:t2_waits_t1_end)
            @record(:t2_waits_t3_begin)
            wait(t3)  # waiting a s1 task in a s2 task
            @record(:t2_waits_t3_end)
        end
    end
    @test fetch_finished(t1) == 111
    @test fetch_finished(t2) == 222
    @test fetch_finished(t3) == 333
end

function test_combination2()
    @testset "$(p.label)" for p in sets2()
        check_combination2(p.s1, p.s2)
    end
end

end  # module


text:
using Distributions
using Gadfly
using Colors

"""
Bayes changes of variables
"""
function drawBayesChangeOfVar()
    xseq = collect(0:0.1:10)

    ginv(x) = 1./(1+exp(-x+5))

    μ = 6
    Σ = 1
    N = 10^6

    rn = Normal(μ, Σ)
    x  = rand(rn, N)

    yheight, ypoints = hist(ginv(x), 51)

    titles=["g"; "p_X"; "p_Y"]
    colors=["blue"; "red"; "green"]

    g   = layer(x=xseq, y=ginv(xseq), Geom.line, Theme(default_color=parse(Colorant, colors[1]), line_width=3px))
    p_x = layer(x=x, Geom.histogram(bincount=50, density=true), Theme(default_color=parse(Colorant, colors[2])))
    p_y = layer(y=ginv(x), Geom.histogram(bincount=50, density=true, orientation=:horizontal),
               Theme(default_color=parse(Colorant, colors[3])))
    hl  = layer(x=[0; μ], y=[ginv(μ); ginv(μ)], Geom.line, Theme(line_width=3px))
    vl  = layer(x=[μ; μ], y=[0, ginv(μ)], Geom.line, Theme(line_width=3px))

    pp = plot(hl, vl, g, p_x, p_y, Guide.manual_color_key("Bayes changes of vars", titles, colors))
end


text:
# ---
# title: 284. Peeking Iterator
# id: problem284
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Design
# link: <https://leetcode.com/problems/peeking-iterator/description/>
# hidden: true
# ---
#
# Given an Iterator class interface with methods: `next()` and `hasNext()`,
# design and implement a PeekingIterator that support the `peek()` operation --
# it essentially peek() at the element that will be returned by the next call to
# next().
#
# **Example:**
#
#
#
#     Assume that the iterator is initialized to the beginning of the list: **[1,2,3]**.
#
#     Call **next()** gets you **1** , the first element in the list.
#     Now you call **peek()** and it returns **2** , the next element. Calling **next()** after that _**still**_ return **2**.
#     You call **next()** the final time and it returns **3** , the last element.
#     Calling **hasNext()** after that should return **false**.
#
#
# **Follow up** : How would you extend your design to be generic and work with
# all types, not just integer?
#
#
## @lc code=start
using LeetCode

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


text:
# Julia wrapper for header: /usr/include/scip/scip_general.h
# Automatically generated using Clang.jl wrap_c

function SCIPversion()
    ccall((:SCIPversion, libscip), Cdouble, ())
end

function SCIPmajorVersion()
    ccall((:SCIPmajorVersion, libscip), Cint, ())
end

function SCIPminorVersion()
    ccall((:SCIPminorVersion, libscip), Cint, ())
end

function SCIPtechVersion()
    ccall((:SCIPtechVersion, libscip), Cint, ())
end

function SCIPsubversion()
    ccall((:SCIPsubversion, libscip), Cint, ())
end

function SCIPprintVersion(scip, file)
    ccall((:SCIPprintVersion, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintBuildOptions(scip, file)
    ccall((:SCIPprintBuildOptions, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintError(retcode)
    ccall((:SCIPprintError, libscip), Cvoid, (SCIP_RETCODE,), retcode)
end

function SCIPcreate(scip)
    ccall((:SCIPcreate, libscip), SCIP_RETCODE, (Ptr{Ptr{SCIP_}},), scip)
end

function SCIPfree(scip)
    ccall((:SCIPfree, libscip), SCIP_RETCODE, (Ptr{Ptr{SCIP_}},), scip)
end

function SCIPgetStage(scip)
    ccall((:SCIPgetStage, libscip), SCIP_STAGE, (Ptr{SCIP_},), scip)
end

function SCIPprintStage(scip, file)
    ccall((:SCIPprintStage, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPgetStatus(scip)
    ccall((:SCIPgetStatus, libscip), SCIP_STATUS, (Ptr{SCIP_},), scip)
end

function SCIPprintStatus(scip, file)
    ccall((:SCIPprintStatus, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPisTransformed(scip)
    ccall((:SCIPisTransformed, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPisExactSolve(scip)
    ccall((:SCIPisExactSolve, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPisPresolveFinished(scip)
    ccall((:SCIPisPresolveFinished, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPhasPerformedPresolve(scip)
    ccall((:SCIPhasPerformedPresolve, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPpressedCtrlC(scip)
    ccall((:SCIPpressedCtrlC, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPisStopped(scip)
    ccall((:SCIPisStopped, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPincludeExternalCodeInformation(scip, name, description)
    ccall((:SCIPincludeExternalCodeInformation, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring), scip, name, description)
end

function SCIPgetExternalCodeNames(scip)
    ccall((:SCIPgetExternalCodeNames, libscip), Ptr{Cstring}, (Ptr{SCIP_},), scip)
end

function SCIPgetExternalCodeDescriptions(scip)
    ccall((:SCIPgetExternalCodeDescriptions, libscip), Ptr{Cstring}, (Ptr{SCIP_},), scip)
end

function SCIPgetNExternalCodes(scip)
    ccall((:SCIPgetNExternalCodes, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPprintExternalCodes(scip, file)
    ccall((:SCIPprintExternalCodes, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end


text:
"A lens into the soul of your program"
module Lens

using Cassette
using Spec

"""
Create a named lens

```julia
```
"""
function lens end

"Apply f to args in context of lens"
function lensapply end

"Evaluate an expression with a lensmap"
function leval end

# Cassette based implementation
# include("cassette.jl")

# TODO: Normal Julia implementation
# include("julia.jl")

include("global.jl")

"""Function call with lens (shorthand for)

```julia
f(x) = lens(:mylens, 2x + 1)
lmap = (mylens = println,)
@lenscall lmap f(31)
```
"""
macro lenscall(lmap, fcall)
  fcall.head == :call || error("Must be function application")
  haskwargs(expr) = expr.args[2] isa Expr && expr.args[2].head == :parameters
  r = if haskwargs(fcall)
    Expr(:call, :lenscall, fcall.args[2], lmap, fcall.args[1], fcall.args[3:end]...)
  else
    Expr(:call, :lenscall, lmap, fcall.args...)
  end
  r
end
@spec call.head == :call "Must be function application"

export lens, lenscall, @lenscall, @leval

end


text:
using ArrayInterface, Test

@test ArrayInterface.ismutable(rand(3))

using StaticArrays
ArrayInterface.ismutable(@SVector [1,2,3]) == false
ArrayInterface.ismutable(@MVector [1,2,3]) == true


text:
using Documenter, DocumenterCitations
using ProximalAlgorithms, ProximalCore
using Literate

bib = CitationBibliography(joinpath(@__DIR__, "references.bib"))

src_path = joinpath(@__DIR__, "src/")

literate_directories = joinpath.(
    src_path,
    [
        "guide",
        "examples",
    ]
)

for directory in literate_directories
    jl_files = filter(p -> endswith(p, ".jl"), readdir(directory; join=true))
    for src in jl_files
        Literate.markdown(src, directory, documenter=true)
    end
end

makedocs(
    bib,
    modules=[ProximalAlgorithms, ProximalCore],
    sitename="ProximalAlgorithms.jl",
    pages=[
        "Home" => "index.md",
        "User guide" => [
            joinpath("guide", "getting_started.md"),
            joinpath("guide", "implemented_algorithms.md"),
            joinpath("guide", "custom_objectives.md"),
            joinpath("guide", "custom_algorithms.md"),
        ],
        "Examples" => [
            joinpath("examples", "sparse_linear_regression.md"),
        ],
        "Bibliography" => "bibliography.md",
    ],
)

deploydocs(
    repo="github.com/JuliaFirstOrder/ProximalAlgorithms.jl.git",
)


text:
###################################
## Simple VI & GS for 2D Gaussian

using PyPlot
using Distributions

function calc_KL(mu1, lambda1, mu2, lambda2)
    D = size(mu1, 1)
    px_lnqx = 0.5 * logdet(lambda2) - 0.5 * ((mu1 - mu2)' * lambda2 * (mu1 - mu2) + trace(lambda2 * inv(lambda1)))
    px_lnpx = 0.5 * logdet(lambda1) - 0.5 * D
    KL = - (px_lnqx - px_lnpx)
    return KL[1]
end

function plot_results(result, truth)
    N = size(result, 1)
    H = Int(ceil(sqrt(N)))
    W = Int(ceil(N / H))
    for i in 1 : H
        for j in 1 : W
            n = (i - 1) * W + j
            if n <= N
                subplot(H, W, n)
                title("$n of $N")
                plot_gaussian(truth[1], truth[2], "b", "\$p(z)\$")
                plot_gaussian(result[n][1], result[n][2], "r", "\$p(z)\$")
            end
        end
    end
end

function plot_lines(X)
    D, N = size(X)
    X_d = zeros(D, 2*N + 1)
    X_d[:,1] = X[:,1]
    for i in 1 : N
        X_d[1, 2*i - 1] = X[1, i]
        X_d[1, 2*i] = X[1, i]
        X_d[2, 2*i] = X[2, i]
        X_d[2, 2*i + 1] = X[2, i]
    end
    plot(X[1,:], X[2,:], "oy")
    plot(X_d[1,1:2*N], X_d[2,1:2*N], "--y")
end

function plot_gaussian(Mu, Sigma, col, label)
    res = 100
    plot(Mu[1], Mu[2], "x", color=col)

    F = eigfact(Sigma)
    vec = F.vectors
    val = F.values
    dw = 2*pi/res
    w = dw * (0 : res)

    c = 1.0
    a = sqrt(c*val[1])
    b = sqrt(c*val[2])
    P1 = a*cos.(w)
    P2 = b*sin.(w)
    P = Mu .+ vec'*vcat(P1', P2')
    plot(P[1, :], P[2, :], "-", color=col, label=label)
end

"""
Variational inference for 2D Gauss.
"""
function main_VI()
    ## creat truth distribution
    D = 2 # dimension
    theta = 2.0*pi/12 # tilt
    A = reshape([cos.(theta), -sin.(theta),
                 sin.(theta), cos.(theta)],
                2, 2)
    mu = [0.0, 0.0]
    lambda = inv(A * inv(reshape([1,0,0,10], 2, 2)) * A')

    ## initialize
    #mu_h = randn(D)
    mu_h = [-0.5, 0.3]
    lambda_h = zeros(D,D)

    ## main iteration
    max_iter = 10
    KL = NaN * Array{Float64, 1}(max_iter)
    result = Array{Any, 1}(max_iter)
    for i in 1 : max_iter
        ## update
        mu_h[1] = mu[1] - inv(lambda[1,1])*lambda[1,2] * (mu_h[2] - mu[2])

        lambda_h[1,1] = lambda[1,1]
        mu_h[2] = mu[2] - inv(lambda[2,2])*lambda[2,1] * (mu_h[1] - mu[1])
        lambda_h[2,2] = lambda[2,2]

        ## calculate KL divergeince
        KL[i] = calc_KL(mu_h, lambda_h, mu, lambda)

        ## store the results
        result[i] = [deepcopy(mu_h), deepcopy(inv(lambda_h))]
    end

    ## visualize results
    figure("result per iteration (VI)")
    clf()
    plot_results(result, (mu, inv(lambda)))

    figure("result (VI)")
    clf()
    plot_gaussian(mu, inv(lambda), "b", "\$p(\\bf{z})\$")
    plot_gaussian(result[end][1], result[end][2], "r", "\$q(\\bf{z})\$")
    xlabel("\$z_1\$", fontsize=20)
    ylabel("\$z_2\$", fontsize=20)
    legend(fontsize=16)

    figure("KL divergence (VI)")
    clf()
    plot(1:max_iter, KL)
    ylabel("KL divergence", fontsize=16)
    xlabel("iteration", fontsize=16)
    show()
end

"""
Gibbs sampling for 2D Gauss.
"""
function main_GS()
    ## creat truth distribution
    D = 2 # dimension
    theta = 2.0*pi/12 # tilt
    A = reshape([cos.(theta), -sin.(theta),
                 sin.(theta), cos.(theta)],
                2, 2)
    mu = [0.0, 0.0]
    #lambda = inv(A * inv(reshape([1,0,0,10], 2, 2)) * A')
    lambda = inv(A * inv(reshape([1,0,0,100], 2, 2)) * A')

    ## initialize
    #max_iter = 1000
    max_iter = 50
    X = randn(D, max_iter)
    mu_h = randn(D)

    ## main iteration
    KL = NaN * Array{Float64, 1}(max_iter)
    for i in 2 : max_iter
        ## update
        mu_h[1] = mu[1] - inv(lambda[1,1])*lambda[1,2] * (X[2,i-1] - mu[2])
        X[1, i] = rand(Normal(mu_h[1], sqrt(inv(lambda[1,1]))))

        mu_h[2] = mu[2] - inv(lambda[2,2])*lambda[2,1] * (X[1,i] - mu[1])
        X[2, i] = rand(Normal(mu_h[2], sqrt(inv(lambda[2,2]))))

        if i > D
            KL[i] = calc_KL(mean(X[:,1:i], 2), inv(cov(X[:,1:i], 2)), mu, lambda)
        end
    end

    ## visualize results
    expt_mu = mean(X, 2)
    expt_Sigma = cov(X, 2)

    figure("samples (GS)")
    clf()
    plot_lines(X)
    plot_gaussian(mu, inv(lambda), "b", "\$p(\\bf{z})\$")
    plot_gaussian(expt_mu, expt_Sigma, "r", "\$q(\\bf{z})\$")
    xlabel("\$z_1\$", fontsize=20)
    ylabel("\$z_2\$", fontsize=20)
    legend(fontsize=16)

    figure("KL divergence (GS)")
    clf()
    plot(1:max_iter, KL)
    ylabel("KL divergence", fontsize=16)
    xlabel("sample size", fontsize=16)
    show()
end

main_VI()
main_GS()

max_stars_repo_name: sammy-suyama/BayesBook

max_issues_repo_name: sammy-suyama/BayesBook

max_forks_repo_name: sammy-suyama/BayesBook


text:
#################### File I/O ####################

function Base.read(name::AbstractString, ::Type{T}) where {T<:Chains}
  c = open(deserialize, name, "r")
  isa(c, T) || throw(TypeError(:open, "read(\"$name\", $T)", T, c))
  return c
end

function Base.write(name::AbstractString, c::Chains)
  open(file -> serialize(file, c), name, "w")
end

function readcoda(output::AbstractString, index::AbstractString)
  out = readdlm(output, Any)
  ind = readdlm(index, Any)

  firstind = ind[:, 2]
  firstiter = out[firstind, 1]
  lastind = ind[:, 3]
  lastiter = out[lastind, 1]

  thin = Int((lastiter[1] - firstiter[1]) / (lastind[1] - firstind[1]))
  window = maximum(firstiter):thin:minimum(lastiter)
  startind = firstind + (first(window) - firstiter) / step(window)
  stopind = lastind - (lastiter - last(window)) / step(window)

  names = AbstractString[ind[:, 1]...]

  value = Array{Float64}(length(window), length(names))
  for i in 1:size(value, 2)
    inds = Int(startind[i]):Int(stopind[i])
    value[:, i] = out[inds, 2]
  end

  Chains(value, start=first(window), thin=step(window), names=names)
end


text:
@testset "Test Problems" begin
    @testset "Unconstrained" begin
        @testset "camel6" begin
            p = DSProblem(2; poll=OrthoMADS(2), objective=DS.camel6, initial_point=[5.0, 5.0])
            Optimize!(p)
            @test isapprox(abs(p.x[1]), 0.0898, atol=1e-4)
            @test isapprox(abs(p.x[2]), 0.7126, atol=1e-4)
            @test isapprox(p.x_cost, -1.0316, atol=1e-4)
            @test p.iteration < 120
        end
        @testset "audet_dennis_benchmark_1" begin
            p = DSProblem(2; poll=OrthoMADS(2), objective=DS.bm_1, initial_point=[-2.1, 1.7])
            Optimize!(p)
            @test isapprox(p.x[1], 0, atol=1e-8)
            @test isapprox(p.x[2], 0, atol=1e-8)
            @test isapprox(p.x_cost, 0, atol=1e-8)
            @test p.iteration < 120
        end
    end
    @testset "Extreme Constraints" begin
        #@testset "camel6" begin
        #    p = DSProblem(2; poll=OrthoMADS(2), objective=DS.camel6, initial_point=[5.0, 5.0])
        #    AddExtremeConstraint(p, x-> x[1] > 0)
        #    Optimize!(p)
        #    @test isapprox(p.x[1], 0.0898, atol=1e-4)
        #    @test isapprox(p.x[2], -0.7126, atol=1e-4)
        #    @test isapprox(p.x_cost, -1.0316, atol=1e-4)
        #    @test p.iteration < 120

        #    p = DSProblem(2; poll=OrthoMADS(2), objective=DS.camel6, initial_point=[-5.0, -5.0])
        #    AddExtremeConstraint(p, x->x[1] < 0)
        #    Optimize!(p)
        #    @test isapprox(p.x[1], -0.0898, atol=1e-4)
        #    @test isapprox(p.x[2], 0.7126, atol=1e-4)
        #    @test isapprox(p.x_cost, -1.0316, atol=1e-4)
        #    @test p.iteration < 120
        #end
    end
    @testset "Progressive Constraints" begin
        #Progressive barrier constraints struggle to converge with LTMADS, therefore use OrthoMADS
        @testset "camel6" begin
            p = DSProblem(2; poll=OrthoMADS(2), objective=DS.camel6, initial_point=[5.0, 5.0])
            AddProgressiveConstraint(p, x -> x[1])
            Optimize!(p)
            @test isapprox(p.x[1], -0.0898, atol=1e-4)
            @test isapprox(p.x[2], 0.7126, atol=1e-4)
            @test isapprox(p.x_cost, -1.0316, atol=1e-4)
            @test p.iteration < 200

            p = DSProblem(2; poll=OrthoMADS(2), objective=DS.camel6, initial_point=[-5.0, -5.0])
            AddProgressiveConstraint(p, x->-x[1])
            Optimize!(p)
            @test isapprox(p.x[1], 0.0898, atol=1e-4)
            @test isapprox(p.x[2], -0.7126, atol=1e-4)
            @test isapprox(p.x_cost, -1.0316, atol=1e-4)
            @test p.iteration < 200
        end
    end
end


text:
mutable struct DotAcc{T}
    s :: T
end

dotAcc(T) = DotAcc{T}(vzero(T))

function add!(acc::DotAcc, x, y)
    SIMDops.@explicit
    acc.s += x * y
end

function add!(acc::DotAcc{T}, x::DotAcc{T}) where {T}
    SIMDops.@explicit
    acc.s += x.s
end

Base.sum(acc::DotAcc{T}) where {T<:Vec}  = DotAcc(vsum(acc.s))
Base.sum(acc::DotAcc{T}) where {T<:Real} = acc.s


text:
test_env = mktempdir()

function with_test_env()
    hooks_path = joinpath(@__DIR__, "..") |> normpath

    """
    begin
        import Pkg
        Pkg.activate("$test_env")
        Pkg.develop(path="$hooks_path")
    end
    """
end

function noerror(cell)
    errored = cell.errored
    if errored
        @show cell.output
    end
    !errored
end

function setcode(cell, code)
    cell.code = code
end


text:
module Tube

function get_fn(argd::Dict)
    radius = float(argd["radius"])
    k = float(argd["k"])
    center = [float(x) for x in argd["center"]]
    power = argd["exponent"]
    function nomalize{T<:AbstractFloat}(xs::Array{T,1})
        if norm(xs) < 1e-7
            return xs
        else
            return xs / norm(xs)
        end
    end
    direction = normalize([float(x) for x in argd["direction"]])

    function distance{T<:AbstractFloat}(point::Array{T,1})
        p = point - center
        return norm(p - dot(p, direction)*direction)
    end
    function fn_compute(coordinates)
        d = distance(coordinates)
        if d <= radius
            return 0.0
        else
            return 0.5 * k * d^power
        end
    end

    return fn_compute
end

export get_fn
end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
using Compat

if @compat is_unix()
  cd("/tmp")
end

info(" running tests in: $(pwd())")

info("starting test arc-twopoints" );               include("arc-twopoints.jl")
info("starting test arrow-arc-test" );              include("arrow-arc-test.jl")
info("starting test arrow-line-test" );             include("arrow-line-test.jl")
info("starting test axes-test" );                   include("axes-test.jl")
info("starting test circletests" );                 include("circletests.jl")
info("starting test clipping-test" );               include("clipping-test.jl")

# this test fails when followed by other tests in the same session, but not on its own.
# Presumably there's some pointer or something which isn't released properly
# info("starting test get_path" ); include("get_path.jl")

info("starting test ellipse test" );               include("ellipse-test.jl")
info("starting test heart-julia" );                include("heart-julia.jl")
info("starting test images_with_alpha" );          include("images-with-alpha.jl")
info("starting test julia_logo_draw_eps" );        include("julia-logo-draw-eps.jl")
info("starting test julia_logo_draw" );            include("julia-logo-draw.jl")
info("starting test luxor-test1") ;                include("luxor-test1.jl")
info("starting test matrix-tests" );               include("matrix-tests.jl")
info("starting test palette_test" );               include("palette_test.jl")
info("starting test pie-test" );                   include("pie-test.jl")
info("starting test point-inside-polygon" );       include("point-inside-polygon.jl")
info("starting test point-intersection" );         include("point-intersection.jl")
info("starting test polygon-centroid-sort-test" ); include("polygon-centroid-sort-test.jl")
info("starting test polygon-test" );               include("polygon-test.jl")
info("starting test polysplit" );                  include("polysplit.jl")
info("starting test pretty-poly-test" );           include("pretty-poly-test.jl")
info("starting test randomsinecurves" );           include("randomsinecurves.jl")
info("starting test sector-test" );                include("sector-test.jl")
info("starting test sierpinski-svg" );             include("sierpinski-svg.jl")
info("starting test sierpinski" );                 include("sierpinski.jl")
info("starting test simplify-polygons" );          include("simplify-polygons.jl")
info("starting test star-test-1" );                include("star-test-1.jl")
info("starting test star-test" );                  include("star-test.jl")
info("starting test test-holes" );                 include("test-holes.jl")
info("starting test test-image" );                 include("test-image.jl")
info("starting test text-alignment" );             include("text-alignment.jl")
info("starting test text-path-clipping" );         include("text-path-clipping.jl")
info("starting test tiling-images" );              include("tiling-images.jl")
info("starting test turtle" );                     include("turtle.jl")

info("all tests finished. Images were saved to $(pwd()).")


text:
using Test
using RegionTrees
using StaticArrays
using LinearAlgebra

include("hyperrectangle.jl")
include("cell.jl")
include("twosarray.jl")
include("asdfs.jl")


text:
include("4.1 Sine/.jl")
include("4.3 Tangent/.jl")
include("4.5 Secant/.jl")
include("4.7 Miscellaneous/.jl")

max_stars_repo_path: src/rules/4 Trig functions/.jl

max_issues_repo_path: src/rules/4 Trig functions/.jl

max_forks_repo_path: src/rules/4 Trig functions/.jl