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

using Markdown
using InteractiveUtils

# ╔═╡ 1bcf25d8-71c2-11eb-05ca-6fe7c4e1b24f
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)

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

# ╔═╡ 96e7042c-73be-11eb-3ac5-9f9998b8c8ff
using LinearAlgebra

# ╔═╡ d17897d0-7540-11eb-1d5e-191caf244165
include("./utils.jl")

# ╔═╡ 9dcffb1c-71c1-11eb-2402-39eaa1f1b5e1
md"""
## Week 4: Ridge Regression (Gradient Descent)
"""

# ╔═╡ c39f9348-71c1-11eb-1691-1d28c6b1d916
md"""
In this notebook we will cover estimating multiple regression weights via gradient descent, performing the following:

  - Write a function to compute the derivative of the regression weights with respect to a single feature
  -  Write gradient descent function to compute the regression weights given an initial weight vector, step size, tolerance, and L2 penalty
"""

# ╔═╡ 3e5d907e-71c2-11eb-2eb3-4375fb21b884
md"""
### Load in house sales data
"""

# ╔═╡ 4c736078-71c2-11eb-225b-816745199225
sales = CSV.File("../../ML_UW_Spec/C02/data/kc_house_test_data.csv";
	header=true) |> DataFrame;

# ╔═╡ 35b5d05a-71cc-11eb-15a9-852c5a534fab
md"""
### Convert to Julia Matrix/Vector
"""

# ╔═╡ a2928dc6-71c2-11eb-3724-29dc1448a0f3
md"""
Now we will write a function that will accept a DataFrame, a list of feature names (e.g. [:sqft_living, :bedrooms') and a target feature e.g. (:price) and will return two things:

 - A matrix whose columns are the desired features plus a constant column (this is how we create an 'intercept')
 - An array containing the values of the output

With this in mind, let's write the `get_data` and `predict_output` functions (moved to `utils.jl`) as follows:

```julia
function get_data(df, features, output)
	df[:, :constant] .= 1.0
	features = [:constant, features...]
	X_matrix = convert(Matrix, select(df, features)) # to get a matrix
	y = df[!, output]                                # => to get a vector
	(X_matrix, y)
end

function predict_output(X::Matrix{T}, weights::Vector{T}) where {T <: Real}
    # assume feature_matrix is a matrix containing the features as columns
	# and weights is a corresponding array
    X * weights
end
```

"""

# ╔═╡ 21635b6e-71c8-11eb-2272-edec98c58339
md"""
### Computing the derivative

We are now going to move to computing the derivative of the regression cost function. Recall that the cost function is the sum over the data points of the squared difference between an observed output and a predicted output, plus the L2 penalty term.

$$Cost(w) = Σ [(prediction - output)^2 + l_2 × Σ_{i=1}^{k} w_i^2]$$

Since the derivative of a sum is the sum of the derivatives, we can take the derivative of the first part (the RSS) as we did in the notebook for the unregularized case in Week 2 and add the derivative of the regularization part. <br />
As we saw, the derivative of the RSS with respect to `w[i]` can be written as:

$$2 × Σ_i [error × [feature_i]]$$

The derivative of the regularization term with respect to `w[i]` is:

$$2 × l_2 \times w_i$$

Summing both, we get

$$2 \times Σ_i [error \times [feature_i] + 2 × l_2 * w_i]$$

That is, the derivative for the weight for feature i is the sum (over data points) of 2 times the product of the error and the feature itself, plus $2 × l_2 * w_i$.

**We will not regularize the constant.**  Thus, in the case of the constant, the derivative is just twice the sum of the errors (without the $2 × l_2 * w_0$ term).

Recall that twice the sum of the product of two vectors is just twice the dot product of the two vectors. Therefore the derivative for the weight for feature_i is just two times the dot product between the values of feature_i and the current errors, plus $2 × l_2 * w_i$.

With this in mind complete the following derivative function which computes the derivative of the weight given the value of the feature (over all data points) and the errors (over all data points).

To decide when to we are dealing with the constant (so we don't regularize it) we added the extra parameter to the call `feature_is_constant` which you should set to `true` when computing the derivative of the constant and `false` otherwise.

"""

# ╔═╡ f377f6d4-71c9-11eb-0249-f3d63e87062b
function feature_derivative_ridge(errors::Vector{T}, feature::Vector{T},
		weights, l2; is_constant=false) where {T <: Real}
	deriv_err = 2 * dot(feature, errors) # ≡ 2 * feature' * errors
	deriv_l2p = l2 * weights
    return is_constant ? deriv_err : deriv_err + deriv_l2p
end

# ╔═╡ 59e3f22a-73bd-11eb-1c0d-150d7d462676
begin
	(ex_features, ex_output) = get_data(sales, [:sqft_living], :price)
	typeof(ex_features), typeof(ex_output)
	ex_weights = Float64[1., 10.]
	test_preds = predict_output(ex_features, ex_weights)
	errors = test_preds - ex_output # prediction errors

	@test feature_derivative_ridge(errors, ex_features[:, 2], ex_weights[2], 1;
		is_constant=false) ≈ sum(errors .* ex_features[:, 2]) * 2. + 20.

	@test feature_derivative_ridge(errors, ex_features[:, 1], ex_weights[1], 1;
		is_constant=true) ≈ sum(errors) * 2.
end

# ╔═╡ 239ff6ca-71cc-11eb-341a-df1453d60a15
md"""
### Gradient Descent [GD]

Now we will write a function that performs a gradient descent. The basic premise is simple. Given a starting point we update the current weights by moving in the negative gradient direction. Recall that the gradient is the direction of increase and therefore the negative gradient is the direction of decrease and we're trying to minimize a cost function.

The amount by which we move in the negative gradient direction is called the *step size* denoted by η. We stop when we are 'sufficiently close' to the optimum. Unlike in Week 2, this time we will set a maximum number of iterations and take gradient steps until we reach this maximum number. If no maximum number is supplied, the maximum should be set 100 by default. (Use default parameter values in Julia.)

With this in mind, write the following gradient descent function below using your derivative function above. For each step in the gradient descent, we update the weight for each feature before computing our stopping criteria.
"""

# ╔═╡ 958f81a6-71cc-11eb-13e8-61dddeb014a7
function regression_gradient_descent(f_matrix::Matrix{T},
		output::Vector{T}, init_weights::Vector{T}, η::T, l2_penalty;
		max_iter=100) where {T <: Real}
    weights = copy(init_weights)
	iter = 0

    while iter < max_iter
		iter += 1
		preds = predict_output(f_matrix, weights) # compute the predictions
        errors = preds .- output                  # compute the errors

		## Update the weights
        for ix ∈ 1:length(weights)                # loop over each weight
            ## Recall that feature_matrix[:, ix] is the feature column
		    ## associated with weights[i]
            deriv_ix = feature_derivative_ridge(errors, f_matrix[:, ix],
				weights[ix],
				l2_penalty;
				is_constant=ix == 1)

			## subtract the step size times the derivative from the CURRENT weight
            weights[ix] -= η * deriv_ix
		end
	end
    weights
end

# ╔═╡ 336f9d62-73c0-11eb-1861-3fb709350f2a
md"""
### Visualizing effect of L2 penalty

The L2 penalty gets its name because it causes weights to have small L2 norms than otherwise. Let's see how large weights get penalized. Let us consider a simple model with 1 feature:
"""

# ╔═╡ e8757854-73c9-11eb-2fb9-c10ac86bf9f0
begin
	sales_train, sales_test = train_test_split(sales; shuffled=false);
	first(sales_train, 3)
end

# ╔═╡ e9114fbe-73cd-11eb-052b-15279aaada23
begin
	s_features = [:sqft_living]
	s_output = :price

	## use: sales_train, sales_test
	(s_feature_mtr, s_output_tr) = get_data(sales_train, s_features, s_output)
	(s_feature_mte, s_output_te) = get_data(sales_test, s_features, s_output)

	s_init_weights = Float64[0., 0.]
	η = 1e-12
	max_iter=1000

	s_l2_penalty = 0.0    # no penalty
	h_l2_penalty = 1e11   # high penalty

	size(s_feature_mtr), size(s_output_tr), size(s_feature_mte), size(s_output_te)
end

# ╔═╡ 47de5248-73c0-11eb-2a0c-85982a16d3a3
s_weights_0_penalty = regression_gradient_descent(s_feature_mtr,
		Vector{Float64}(s_output_tr), s_init_weights, η, s_l2_penalty;
		max_iter)

# ╔═╡ 5cb2aee6-73c3-11eb-0eb0-cfbfcb4c145a
s_weights_h_penalty = regression_gradient_descent(s_feature_mtr,
		Vector{Float64}(s_output_tr), s_init_weights, η, h_l2_penalty;
		max_iter)

# ╔═╡ c20a95dc-73ce-11eb-072e-cfdf28bb8412
(s_weights_0_penalty, s_weights_h_penalty)

# ╔═╡ 87b80b2c-73c3-11eb-3175-0dedf19b1858
begin
	scatter(s_feature_mtr, s_output_tr, color=[:orange], marker=".")
	plot!(s_feature_mtr, predict_output(s_feature_mtr, s_weights_0_penalty), color=[:blue])
	plot!(s_feature_mtr, predict_output(s_feature_mtr, s_weights_h_penalty), color=[:green])
end

# ╔═╡ a2b905e0-73cb-11eb-0aef-79f5e6941f30
md"""
Compute the RSS on the *test* data for the following three sets of weights:
  1. The initial weights (all zeros)
  1. The weights learned with no regularization
  1. The weights learned with high regularization

Which weights perform best?
"""

# ╔═╡ b8f7c40e-73cb-11eb-1e61-8f84439c8cb2
md"""

moved to `utils.jl`

```julia
function (X::Matrix{T}, y, weights::Vector{T})
    preds = predict_output(X, weights)
	sum((preds .- y) .^ 2)
end
```
"""

# ╔═╡ cd05b5fe-73cc-11eb-35a4-952fa79374f4
begin
	## RSS initial weigths
	v0 = calc_rss(s_feature_mte, s_output_te, s_init_weights)

	## RSS weigths no regularization
	v1 = calc_rss(s_feature_mte, s_output_te, s_weights_0_penalty)

	## RSS weigths with high regularization
	v2 = calc_rss(s_feature_mte, s_output_te, s_weights_h_penalty)

	with_terminal() do
		@printf("rss init. weights:     %1.5e\n", v0)
		@printf("rss no reg. weights:   %1.5e\n", v1)
		@printf("rss high reg. weights: %1.5e\n", v2)
	end
end

# ╔═╡ 92029980-73cd-11eb-1de5-9fd4c888d352
md"""

**Quiz Questions**

1. What is the value of the coefficient for `sqft_living` that you learned with no regularization, rounded to 1 decimal place?  What about the one with high regularization?
   - 263.0 (no regularization), 68.7 (with high regularization)

2. Comparing the lines you fit with the with no regularization versus high regularization, which one is steeper?
   - the blue line (with no regularization)

3. What are the RSS on the test data for each of the set of weights above (initial, no regularization, high regularization)?
   - cf. above
"""

# ╔═╡ 27319a50-73d1-11eb-13ee-5b868486e887
md"""
### Running a multiple regression with L2 penalty

Let us now consider a model with 2 features: `[:sqft_living, :sqft_living15]`.

"""

# ╔═╡ 4def015a-73d1-11eb-07be-af419657dd3d
begin
	m_features = [:sqft_living, :sqft_living15]
	m_output = :price

	## use: sales_train, sales_test
	(m_feature_mtr, m_output_tr) = get_data(sales_train, m_features, m_output)
	(m_feature_mte, m_output_te) = get_data(sales_test, m_features, m_output)

	m_init_weights = Float64[0., 0., 0.]

	## As above for the following:
	# η = 1e-12
	# max_iter=1000
	# s_l2_penalty = 0.0    # no penalty
	# h_l2_penalty = 1e11   # high penalty
end

# ╔═╡ 505bc8cc-73d1-11eb-104b-f1fe8928f397
m_weights_0_penalty = regression_gradient_descent(m_feature_mtr,
		Vector{Float64}(m_output_tr), m_init_weights, η, s_l2_penalty;
		max_iter)

# ╔═╡ 502fc1c8-73d1-11eb-3102-4303736bcc5f
m_weights_h_penalty = regression_gradient_descent(m_feature_mtr,
		Vector{Float64}(m_output_tr), m_init_weights, η, h_l2_penalty;
		max_iter)

# ╔═╡ 4ff2dc02-73d1-11eb-333d-07031c47234b
begin
	## RSS initial weigths
	mv0 = calc_rss(m_feature_mte, m_output_te, m_init_weights)

	## RSS weigths no regularization
	mv1 = calc_rss(m_feature_mte, m_output_te, m_weights_0_penalty)

	## RSS weigths with high regularization
	mv2 = calc_rss(m_feature_mte, m_output_te, m_weights_h_penalty)

	with_terminal() do
		@printf("rss init. weights:     %1.5e\n", mv0)
		@printf("rss no reg. weights:   %1.5e\n", mv1)
		@printf("rss high reg. weights: %1.5e\n", mv2)
	end
end

# ╔═╡ 4fdf127a-73d1-11eb-1b41-67a0caaa1eca
begin
	## house price for the 1st house in the test set using the no regularization
	pred_0_no_reg = predict_output(m_feature_mte, m_weights_0_penalty)[1]

	## house price for the 1st house in the test set using high regularization
	pred_0_high_reg = predict_output(m_feature_mte, m_weights_h_penalty)[1]

	with_terminal() do
		@printf("no reg. price: %6.2f / price rounded to nearest dollar: %6.2f\n",
			pred_0_no_reg, round(pred_0_no_reg; digits=0))

		@printf("high reg. price: %6.2f / price rounded to nearest dollar: %6.2f\n",
			pred_0_high_reg, round(pred_0_high_reg; digits=0))

		@printf("actual price: %6.2f\n", m_output_te[1])
	end
end

# ╔═╡ 4fb0a88e-73d1-11eb-3eea-339d461b9e17
abs(pred_0_high_reg - m_output_te[1]), abs(pred_0_no_reg - m_output_te[1])

# ╔═╡ 4f8574ac-73d1-11eb-2b52-dd15c43a1c73
md"""
**Quiz Questions**

 1. What is the value of the coefficient for `sqft_living` that you learned with no regularization, rounded to 1 decimal place?  What about the one with high regularization?
   - 189.3 (no reg), 55.2 (high reg)

 2. What are the RSS on the test data for each of the set of weights above (initial, no regularization, high regularization)?
   -  cf. above

 3. We make prediction for the first house in the test set using two sets of weights (no regularization vs high regularization). Which weights make better prediction *for that particular house*?
   - price from no reg. model

"""

# ╔═╡ Cell order:
# ╟─9dcffb1c-71c1-11eb-2402-39eaa1f1b5e1
# ╟─c39f9348-71c1-11eb-1691-1d28c6b1d916
# ╠═1bcf25d8-71c2-11eb-05ca-6fe7c4e1b24f
# ╠═3e5d907e-71c2-11eb-2eb3-4375fb21b884
# ╠═4c736078-71c2-11eb-225b-816745199225
# ╠═d17897d0-7540-11eb-1d5e-191caf244165
# ╟─35b5d05a-71cc-11eb-15a9-852c5a534fab
# ╟─a2928dc6-71c2-11eb-3724-29dc1448a0f3
# ╟─21635b6e-71c8-11eb-2272-edec98c58339
# ╠═96e7042c-73be-11eb-3ac5-9f9998b8c8ff
# ╠═f377f6d4-71c9-11eb-0249-f3d63e87062b
# ╠═59e3f22a-73bd-11eb-1c0d-150d7d462676
# ╟─239ff6ca-71cc-11eb-341a-df1453d60a15
# ╠═958f81a6-71cc-11eb-13e8-61dddeb014a7
# ╟─336f9d62-73c0-11eb-1861-3fb709350f2a
# ╠═e8757854-73c9-11eb-2fb9-c10ac86bf9f0
# ╠═e9114fbe-73cd-11eb-052b-15279aaada23
# ╠═47de5248-73c0-11eb-2a0c-85982a16d3a3
# ╠═5cb2aee6-73c3-11eb-0eb0-cfbfcb4c145a
# ╠═c20a95dc-73ce-11eb-072e-cfdf28bb8412
# ╠═87b80b2c-73c3-11eb-3175-0dedf19b1858
# ╟─a2b905e0-73cb-11eb-0aef-79f5e6941f30
# ╟─b8f7c40e-73cb-11eb-1e61-8f84439c8cb2
# ╠═cd05b5fe-73cc-11eb-35a4-952fa79374f4
# ╟─92029980-73cd-11eb-1de5-9fd4c888d352
# ╟─27319a50-73d1-11eb-13ee-5b868486e887
# ╠═4def015a-73d1-11eb-07be-af419657dd3d
# ╠═505bc8cc-73d1-11eb-104b-f1fe8928f397
# ╠═502fc1c8-73d1-11eb-3102-4303736bcc5f
# ╠═4ff2dc02-73d1-11eb-333d-07031c47234b
# ╠═4fdf127a-73d1-11eb-1b41-67a0caaa1eca
# ╠═4fb0a88e-73d1-11eb-3eea-339d461b9e17
# ╟─4f8574ac-73d1-11eb-2b52-dd15c43a1c73


text:
reshape(parent::AbstractQuasiArray{T,N}, ndims::Val{N}) where {T,N} = parent
function reshape(parent::AbstractQuasiArray, ndims::Val{N}) where N
    reshape(parent, rdims(Val(N), axes(parent)))
end

ext: jl


text:
function update_omega!(s::StateFS, c::ConstantsFS, d::DataFS, t::TunersFS)
  Q = length(s.omega)
  for q in 1:Q
    update_omega!(q, s, c, d, t)
  end
end

function update_omega!(q::Integer, s::StateFS, c::ConstantsFS, d::DataFS,
                       t::TunersFS)

  function loglike(om_q::Float64)::Float64
    om = copy(s.omega)  # make a copy
    om[q] = om_q
    p = compute_p(d.X, om)

    ll = 0.
    for i in 1:d.data.I
      lpdf(rik::Bool)::Float64 = logpdf(Bernoulli(p[i]), rik)
      for k in 1:c.constants.K
        ll += lpdf(s.r[i, k])
      end
    end

    return ll
  end

  logprior(om_q::Float64)::Float64 = logpdf(c.omega_prior, om_q)
  log_prob(om_q::Float64)::Float64 = loglike(om_q) + logprior(om_q)

  s.omega[q] = MCMC.metropolisAdaptive(s.omega[q], log_prob, t.omega[q])
end


text:
# Autogenerated wrapper script for Tk_jll for aarch64-linux-musl
export libtk, wish

using Tcl_jll
using Xorg_libXft_jll
JLLWrappers.@generate_wrapper_header("Tk")
JLLWrappers.@declare_library_product(libtk, "libtk8.6.so")
JLLWrappers.@declare_executable_product(wish)
function __init__()
    JLLWrappers.@generate_init_header(Tcl_jll, Xorg_libXft_jll)
    JLLWrappers.@init_library_product(
        libtk,
        "lib/libtk8.6.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        wish,
        "bin/wish8.6",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module VariableSelector

using DataFrames
using Statistics
using StatsBase
using HypothesisTests
export variable_importance

# Function for cat/con
# The test for such combo is "KruskalWallis"
function gamma(cat_var::Vector{String}, con_var::Vector{<: Number})
    data = DataFrame(levels = cat_var, values = con_var)
    levels = data.levels |> Set |> collect
    vec_list = Vector{Float64}[]
    for i in levels
      filtered_df = data[data.levels .== i, :]
      filtered_vec = filtered_df.values
      push!(vec_list, filtered_vec)
    end
    p = KruskalWallisTest(vec_list...) |> pvalue
    p = round(p, digits  = 3)
    return (p, "KruskalWallis")
  end

  # Function for con/cat
  # The test for such combo is "KruskalWallis"
  function gamma(con_var::Vector{<: Number}, cat_var::Vector{String})
    gamma(cat_var, con_var)
  end

  # Function for cat/cat
  # The test for such combo is "ChiSquared"
  function gamma(cat_var1::Vector{String}, cat_var2::Vector{String})
    data = DataFrame(levels1 = cat_var1, levels2 = cat_var2)
    ul1 = data.levels1 |> Set |> collect |> sort
    ul2 = data.levels2 |> Set |> collect |> sort
    output = zeros(Int64, length(ul1), length(ul2))
    for i in eachindex(ul1), j in eachindex(ul2)
      tv = (data.levels1 .== ul1[i]) .& (data.levels2 .== ul2[j])
      output[i,j] = sum(tv)
    end
    p = ChisqTest(output) |> pvalue
    return (p, "ChiSquared")
  end

  # Function for con/con
  function gamma(con_var1::Vector{<: Number}, con_var2::Vector{<: Number})
    th = 20 # Treshold for if vector should be category or not

    # Get the # of levels for each vector
    l1 = con_var1 |> Set |> length
    l2 = con_var2 |> Set |> length

    # Get the string vector of each vector
    con_var1_st = string.(con_var1)
    con_var2_st = string.(con_var2)

    # Treat vectors with low number of levels as categories and apply the right test
    if l1 <= th && l2 <= th
      gamma(con_var1_st, con_var1_st) # Apply ChiSquared test
    elseif l1 <= th && l2 > th
      gamma(con_var1_st, con_var2)  # Apply "KruskalWallis" test
    elseif l1 > th && l2 <= th
      gamma(con_var1, con_var2_st) # Apply "KruskalWallis" test
    else
      p = CorrelationTest(con_var1, con_var2) |> pvalue # Apply "Correlation" test
      return (p, "Correlation")
    end
  end

  # Finds the most important variables against the target
  function variable_importance(data::DataFrame, target::Symbol)
    X = select(data, Not(target)) # Get DataFrame X of x variables
    y = data[!, target] # Obtain y vector
    col_names = names(X) # Name vector of X dataframe
    output = DataFrame(feature = String[], test = String[], pvalue = Float64[])
    for i in col_names
      x = X[!, i]
      J = gamma(x, y)
      push!(output, [i, J[2], J[1]])
    end
    tv = (output.pvalue .< 0.05)
    results = output[tv, :]
    return results
  end

end

max_stars_repo_name: ronbactawar/VariableSelection

max_issues_repo_name: ronbactawar/VariableSelection

max_forks_repo_name: ronbactawar/VariableSelection


text:
# Structure for VMLM of OptimPack, see OptimPack.jl
# https://github.com/emmt/OptimPack.jl

ls = OptimPack.MoreThuenteLineSearch(ftol = 1e-8, gtol = 0.95)
scl = OptimPack.SCALING_OREN_SPEDICATO
gat = 0
grt = 1e-3
vt = false
memsize = 100
mxvl = 1000
mxtr = 1000
stpmn = 1e-20
stpmx = 1e+20


text:
#///////////////////////////////////////
#// File Name: nelder_mead_bilevel_optimization.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2020/11/05
#// Description: Nelder-Mead Simplex Bilevel Optimization algorithm
#///////////////////////////////////////

using Printf

"""
    NelderMeadBilevelOptimizationSolver(kwargs...)

RAT iLQR++ (i.e. Nelder-Mead Simplex Method + iLEQG) Solver.

# Optional Keyword Arguments
## iLEQG Solver Parameters
- `μ_min_ileqg::Float64` -- minimum value for Hessian regularization parameter `μ` (> 0).
  Default: `1e-6`.
- `Δ_0_ileqg::Float64` -- minimum multiplicative modification factor (> 0) for `μ`.
  Default: `2.0`.
- `λ_ileqg::Float64` -- multiplicative modification factor in (0, 1) for line search
  step size `ϵ`. Default: `0.5`.
- `d_ileqg::Float64` -- convergence error norm threshold (> 0). If the maximum l2
  norm of the change in nominal control over the horizon is less than `d`, the
  solver is considered to be converged. Default: `1e-2`.
- `iter_max_ileqg::Int64` -- maximum iteration number. Default: 100.
- `ϵ_init_ileqg::Float64` -- initial step size in (`ϵ_min`, 1] to start the
  backtracking line search with. If `adaptive_ϵ_init` is `true`, then this
  value is overridden by the solver's adaptive initialization functionality
  after the first iLEQG iteration. If `adaptive_ϵ_init` is `false`, the
  specified value of `ϵ_init` is used across all the iterations as the initial
  step size. Default:`1.0`.
- `adaptive_ϵ_init_ileqg::Bool` -- if `true`, `ϵ_init` is adaptively changed based on
  the last step size `ϵ` of the previous iLEQG iteration. Default: `false`.
   - If the first line search iterate `ϵ_init_prev` in the previous iLEQG
     iteration is successful, then `ϵ_init` for the next iLEQG iteration is set
     to `ϵ_init = ϵ_init_prev / λ` so that the initial line search step increases.
   - Otherwise `ϵ_init = ϵ_last` where `ϵ_last` is the line search step accepted
     in the previous iLEQG iteration.
- `ϵ_min_ileqg::Float64` -- minimum value of step size `ϵ` to terminate the line
  search. When `ϵ_min` is reached, the last candidate nominal trajectory is accepted
  regardless of the Armijo condition and the current iLEQG iteration is
  finished. Default: `1e-6`.
- `f_returns_jacobian::Bool` -- if `true`, Jacobian matrices of the dynamics function
  are user-provided. This can reduce computation time since automatic
  differentiation is not used. Default: `false`.

## Nelder-Mead Simplex Solver Parameters
- `α::Float64` -- reflection parameter. Default: `1.0`.
- `β::Float64` -- expansion parameter. Default: `2.0`.
- `γ::Float64` -- contraction parameter. Default: `0.5`.
- `ϵ::Float64` -- convergence parameter. The algorithm is said to have convergeced
   if the standard deviation of the objective values at the vertices of the simplex
   is below `ϵ`. Default: `1e-2`.
- `λ::Float64` -- multiplicative modification factor in (0, 1) for `θ_high_init` and
  `θ_low_init`, which is repeatedly applied in case the objective value is infinity
  until a feasible region is find. Default: `0.5`.
- `θ_high_init::Float64` -- Initial guess for `θ_high`. Default: `3.0`.
- `θ_low_init::Float64` -- Initial guess for `θ_low`. Default: `1e-8`.
- `iter_max::Int64` -- maximum iteration number. Default: `100`.

# Notes
- The Nelder-Mead Simplex method maintains a 1D simplex (i.e. a line segment that
  consists of 2 points, `θ_high` and `θ_low`) to search for the optimal risk-sensitivity
  parameter `θ`. `θ_high` and `θ_low` refer to the verteces of the simplex with the highest
  and the lowest objective values, respectively.
- The initial guesses `θ_high_init` and `θ_low_init`, which may be modified during optimization,
  are stored internally in the solver and carried over to the next call to `solve!`.
"""
mutable struct NelderMeadBilevelOptimizationSolver
    # ileqg solver parameters
    μ_min_ileqg::Float64           # Minimum damping parameter for regularization
    Δ_0_ileqg::Float64             # Minimum Modification factor for μ
    λ_ileqg::Float64               # Multiplicative factor for line search step parameter in (0, 1)
    d_ileqg::Float64               # Convergence error norm thresholds
    iter_max_ileqg::Int64          # Maximum iteration
    # β_ileqg::Float64               # Armijo search condition
    ϵ_init_auto_ileqg::Bool        # Automatic initialization of ϵ_init from the previous iLEQG iteration.
    ϵ_init_ileqg::Float64          # Initial step size for backtracking line search
    ϵ_min_ileqg::Float64           # Minimum step size for backtracking line search
    f_returns_jacobian::Bool       # Whether the dynamics funcion f also returns jacobians or not

    # Nelder-Mead Simplex solver parameters
    α::Float64           # Reflection parameter
    β::Float64           # Expansion parameter
    γ::Float64           # Contraction parameter
    ϵ::Float64           # Convergence parameter
    λ::Float64           # Multiplicative factor for θ_high and θ_low in (0, 1)
    iter_max::Int64      # Maximum iteration

    # Nelder-Mead Simplex solver mutable parameters
    θ_high_init::Float64
    θ_low_init::Float64
    iter_current::Int64  # Current Nelder-Mead iteration
    θ_high::Float64
    θ_low::Float64
    c_high::Union{Nothing, Float64}
    c_low::Union{Nothing, Float64}
end

function NelderMeadBilevelOptimizationSolver(;μ_min_ileqg=1e-6,
                                              Δ_0_ileqg=2.0,
                                              λ_ileqg=0.5,
                                              d_ileqg=1e-2,
                                              iter_max_ileqg=100,
                                              # β_ileqg=1e-4,
                                              adaptive_ϵ_init_ileqg=false,
                                              ϵ_init_ileqg=1.0,
                                              ϵ_min_ileqg=1e-6,
                                              f_returns_jacobian=false,
                                              α=1.0,
                                              β=2.0,
                                              γ=0.5,
                                              ϵ=1e-2,
                                              λ=0.5,
                                              iter_max=100,
                                              θ_high_init=3.0,
                                              θ_low_init=1e-8)
    θ_high, θ_low = θ_high_init, θ_low_init;
    iter_current = 0

    return NelderMeadBilevelOptimizationSolver(μ_min_ileqg, Δ_0_ileqg, λ_ileqg, d_ileqg,
                                               iter_max_ileqg, adaptive_ϵ_init_ileqg,
                                               ϵ_init_ileqg, ϵ_min_ileqg, f_returns_jacobian,
                                               α, β, γ, ϵ, λ, iter_max, θ_high_init, θ_low_init,
                                               iter_current, θ_high, θ_low, nothing, nothing)
end

"""
Compute Nelder-Mead cost on a worker process
"""
function compute_cost_worker(nm_solver::NelderMeadBilevelOptimizationSolver,
                             problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                             x::Vector{Float64}, u_array::Vector{Vector{Float64}},
                             θ::Float64, kl_bound::Float64)
    ileqg = ILEQGSolver(problem,
                        μ_min=nm_solver.μ_min_ileqg,
                        Δ_0=nm_solver.Δ_0_ileqg,
                        λ=nm_solver.λ_ileqg,
                        d=nm_solver.d_ileqg,
                        iter_max=nm_solver.iter_max_ileqg,
                        # β=nm_solver.β_ileqg,
                        adaptive_ϵ_init=nm_solver.ϵ_init_auto_ileqg,
                        ϵ_init=nm_solver.ϵ_init_ileqg,
                        ϵ_min=nm_solver.ϵ_min_ileqg,
                        f_returns_jacobian=nm_solver.f_returns_jacobian)
    #initialize!(ileqg, problem, x, u_array, θ)
    cost = 0.0
    try
        cost = solve!(ileqg, problem, x, u_array, θ=θ, verbose=false)[4] +
               kl_bound/θ;
    catch
        cost = Inf
    end
    return cost
end

"""
Initialize NelderMeadBilevelOptimization solver
"""
function initialize!(nm_solver::NelderMeadBilevelOptimizationSolver)
    nm_solver.iter_current = 0;
    nm_solver.θ_low = nm_solver.θ_low_init;
    nm_solver.θ_high = nm_solver.θ_high_init;
end

"""
Single iteration of NelderMeadBilevelOptimization solver
"""
function step!(nm_solver::NelderMeadBilevelOptimizationSolver,
               problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
               x::Vector{Float64},
               u_array::Vector{Vector{Float64}},
               kl_bound::Float64,
               verbose=true)
    nm_solver.iter_current += 1;
    if verbose
        println("**Nelder-Mead iteration $(nm_solver.iter_current)")
    end
    if nm_solver.c_high < nm_solver.c_low
        nm_solver.θ_low, nm_solver.θ_high = nm_solver.θ_high, nm_solver.θ_low;
        nm_solver.c_low, nm_solver.c_high = nm_solver.c_high, nm_solver.c_low;
    end
    if verbose
        println("****(θ_low, c_low) == ($(round(nm_solver.θ_low, digits=4)), $(round(nm_solver.c_low, digits=3)))")
        println("****(θ_high, c_high) == ($(round(nm_solver.θ_high, digits=4)), $(round(nm_solver.c_high, digits=3)))")
    end

    θ_m = nm_solver.θ_low;
    # reflection
    θ_r = θ_m + nm_solver.α*(θ_m - nm_solver.θ_high);
    θ_r = max(nm_solver.θ_low_init, θ_r)
    c_r = compute_cost_worker(nm_solver, problem, x, u_array, θ_r, kl_bound);
    if verbose
        println("****Reflection point: (θ_r, c_r) == ($(round(θ_r, digits=4)), $(round(c_r, digits=3)))")
    end

    if c_r < nm_solver.c_low
        # expansion
        θ_e = θ_m + nm_solver.β*(θ_r - θ_m);
        θ_e = max(nm_solver.θ_low_init, θ_e)
        c_e = compute_cost_worker(nm_solver, problem, x, u_array, θ_e, kl_bound);
        if verbose
            println("****Expansion point: (θ_e, c_e) == ($(round(θ_e, digits=4)), $(round(c_e, digits=3)))")
        end
        if c_e < c_r
            if verbose
                println("****(θ_high, c_high) <-- (θ_e, c_e)")
            end
            nm_solver.θ_high = θ_e;
            nm_solver.c_high = c_e;
        else
            if verbose
                println("****(θ_high, c_high) <-- (θ_r, c_r)")
            end
            nm_solver.θ_high = θ_r;
            nm_solver.c_high = c_r;
        end
    else
        if verbose
            println("****(θ_high, c_high) <-- (θ_r, c_r)")
        end
        if c_r < nm_solver.c_high
            nm_solver.θ_high = θ_r;
            nm_solver.c_high = c_r;
        end
        # contraction
        θ_c = θ_m + nm_solver.γ*(nm_solver.θ_high - θ_m);
        θ_c = max(nm_solver.θ_low_init, θ_c)
        c_c = compute_cost_worker(nm_solver, problem, x, u_array, θ_c, kl_bound);
        if verbose
            println("****Contraction point: (θ_c, c_c) == ($(round(θ_c, digits=4)), $(round(c_c, digits=3)))")
        end
        if c_c > nm_solver.c_high
            nm_solver.θ_high = (nm_solver.θ_high + nm_solver.θ_low)/2;
            nm_solver.c_high = compute_cost_worker(nm_solver, problem, x, u_array, nm_solver.θ_high, kl_bound);
            if verbose
                println("****Shrinking θ_high: (θ_high, c_high) == ($(round(nm_solver.θ_high, digits=4)), $(round(nm_solver.c_high, digits=3)))")
            end
        else
            if verbose
                println("****(θ_high, c_high) <-- (θ_c, c_c)")
            end
            nm_solver.θ_high = θ_c;
            nm_solver.c_high = c_c;
        end
    end
end

"""
    solve!(nm_solver::NelderMeadBilevelOptimizationSolver,
    problem::FiniteHorizonRiskSensitiveOptimalControlProblem, x_0::Vector{Float64},
    u_array::Vector{Vector{Float64}}; kl_bound::Float64, verbose=true)

Given `problem` and `nm_solver` (i.e. a RAT iLQR++ Solver), solve distributionally robust
control with current state `x_0` and nominal control schedule `u_array = [u_0, ..., u_{N-1}]`
under the KL divergence bound of `kl_bound` (>= 0).

# Return Values (Ordered)
- `θ_opt::Float64` -- optimal risk-sensitivity parameter.
- `x_array::Vector{Vector{Float64}}` -- nominal state trajectory `[x_0,...,x_N]`.
- `l_array::Vector{Vector{Float64}}` -- nominal control schedule `[l_0,...,l_{N-1}]`.
- `L_array::Vector{Matrix{Float64}}` -- feedback gain schedule `[L_0,...,L_{N-1}]`.
- `value::Float64` -- optimal cost-to-go (i.e. objective value) found by the solver.

# Notes
- Returns a time-varying affine state-feedback policy `π_k` of the form
  `π_k(x) = L_k(x - x_k) + l_k`.
- If `kl_bound` is 0.0, the solver reduces to iLQG.
"""
function solve!(nm_solver::NelderMeadBilevelOptimizationSolver,
                problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                x_0::Vector{Float64}, u_array::Vector{Vector{Float64}};
                kl_bound::Float64, verbose=true)
    @assert kl_bound >= 0 "KL Divergence Bound must be non-negative"
    initialize!(nm_solver);
    if kl_bound > 0
        if isnothing(nm_solver.c_high)
            while true
                nm_solver.c_high = compute_cost_worker(nm_solver, problem, x_0, u_array, nm_solver.θ_high, kl_bound)
                if !isinf(nm_solver.c_high)
                    break;
                else
                    nm_solver.θ_high *= nm_solver.λ;
                    nm_solver.θ_high_init *= nm_solver.λ;
                end
            end
        end
        if isnothing(nm_solver.c_low)
            while true
                nm_solver.c_low = compute_cost_worker(nm_solver, problem, x_0, u_array, nm_solver.θ_low, kl_bound)
                if !isinf(nm_solver.c_low)
                    break;
                else
                    nm_solver.θ_low *= nm_solver.λ;
                    nm_solver.θ_low_init *= nm_solver.λ;
                end
            end
        end

        while true
            step!(nm_solver, problem, x_0, u_array, kl_bound, verbose)
            # compute stdev
            c_mean = (nm_solver.c_low + nm_solver.c_high)/2;
            stdev = sqrt(0.5*((nm_solver.c_high - c_mean)^2 + (nm_solver.c_low - c_mean)^2));
            stdev_str = @sprintf "%3.3f" stdev
            if stdev < nm_solver.ϵ
                if verbose
                    println("Nelder-Mead Converged. stdev == $(stdev_str)")
                end
                break;
            end
            if nm_solver.iter_current == nm_solver.iter_max
                if verbose
                    println("Maximum iteration number reached. stdev == $(stdev_str)")
                end
                break;
            end
        end
        θ_opt = nm_solver.θ_low;
        if verbose
            println("θ_opt == $(round(θ_opt, digits=4))")
        end
        #nm_solver.θ_high_init = nm_solver.θ_high;
    else
        # kl_bound is 0. In this case the optimal controller is the risk-neutral (iLQG) controller.
        θ_opt = 0.0;
    end
    ileqg_solver = ILEQGSolver(problem,
                               μ_min=nm_solver.μ_min_ileqg,
                               Δ_0=nm_solver.Δ_0_ileqg,
                               λ=nm_solver.λ_ileqg,
                               d=nm_solver.d_ileqg,
                               iter_max=nm_solver.iter_max_ileqg,
                               # β=nm_solver.β_ileqg,
                               adaptive_ϵ_init=nm_solver.ϵ_init_auto_ileqg,
                               ϵ_init=nm_solver.ϵ_init_ileqg,
                               ϵ_min=nm_solver.ϵ_min_ileqg,
                               f_returns_jacobian=nm_solver.f_returns_jacobian)
    #initialize!(ileqg_solver, problem, x_0, u_array, θ_opt)
    x_array, l_array, L_array, value, ~ = solve!(ileqg_solver, problem, x_0, u_array, θ=θ_opt, verbose=verbose)
    if kl_bound > 0
        return θ_opt, x_array, l_array, L_array, value + kl_bound/θ_opt
    else
        return θ_opt, x_array, l_array, L_array, value
    end
end


text:
using ITensors.LazyApply: Applied, AppliedTupleVector

function rrule(f::Type{<:Applied}, x1, x2::Tuple)
  y = f(x1, x2)
  function Applied_pullback(ȳ)
    x̄1 = NoTangent()
    x̄2 = ȳ.args
    return (NoTangent(), x̄1, x̄2)
  end
  return (y, Applied_pullback)
end

function rrule(f::Type{<:AppliedTupleVector}, x1::Vector)
  y = f(x1)
  function Applied_pullback(ȳ)
    x̄1 = ȳ.args[1]
    return (NoTangent(), x̄1)
  end
  function Applied_pullback(ȳ::Vector)
    x̄1 = ȳ
    return (NoTangent(), x̄1)
  end
  function Applied_pullback(ȳ::ZeroTangent)
    x̄1 = ȳ
    return (NoTangent(), x̄1)
  end
  return (y, Applied_pullback)
end


text:
using Fed
using Test

# make sure that everything can start up correctly
include("startup/server.jl")
include("startup/client.jl")

# test the serialization/deserialization
include("payload/serde.jl")


text:
using LazyArrays, ArrayLayouts, LinearAlgebra, FillArrays, StaticArrays, Test
import LazyArrays: BroadcastLayout, arguments, LazyArrayStyle
import Base: broadcasted

@testset "Broadcasting" begin
    @testset "BroadcastArray" begin
        a = randn(6)
        b = BroadcastArray(exp, a)
        @test BroadcastArray(b) == BroadcastVector(b) == b

        @test b ==  Vector(b) == exp.(a)
        @test b[2:5] isa BroadcastVector
        @test b[2:5] == exp.(a[2:5])

        @test exp.(b) isa BroadcastVector
        @test b .+ SVector(1,2,3,4,5,6) isa BroadcastVector
        @test SVector(1,2,3,4,5,6) .+ b isa BroadcastVector

        A = randn(6,6)
        B = BroadcastArray(exp, A)

        @test Matrix(B) == exp.(A)
        @test B[1] == exp(A[1,1])
        @test B[7] == exp(A[1,2])

        C = BroadcastArray(+, A, 2)
        @test C == A .+ 2
        D = BroadcastArray(+, A, C)
        @test D == A + C

        @test sum(B) ≈ sum(exp, A)
        @test sum(C) ≈ sum(A .+ 2)
        @test prod(B) ≈ prod(exp, A)
        @test prod(C) ≈ prod(A .+ 2)

        x = Vcat([3,4], [1,1,1,1,1], 1:3)
        @test x .+ (1:10) isa Vcat
        @test (1:10) .+ x isa Vcat
        @test x + (1:10) isa Vcat
        @test (1:10) + x isa Vcat
        @test x .+ (1:10) == (1:10) .+ x == (1:10) + x == x + (1:10) == Vector(x) + (1:10)

        @test exp.(x) isa Vcat
        @test exp.(x) == exp.(Vector(x))
        @test x .+ 2 isa Vcat
        @test (x .+ 2).args[end] ≡ x.args[end] .+ 2 ≡ 3:5
        @test x .* 2 isa Vcat
        @test 2 .+ x isa Vcat
        @test 2 .* x isa Vcat

        A = Vcat([[1 2; 3 4]], [[4 5; 6 7]])
        @test A .+ Ref(I) == Ref(I) .+ A == Vcat([[2 2; 3 5]], [[5 5; 6 8]])

        @test BroadcastArray(*,1.1,[1 2])[1] == 1.1

        B = BroadcastArray(*, Diagonal(randn(5)), randn(5,5))
        @test B == broadcast(*,B.args...)
        @test colsupport(B,1) == rowsupport(B,1) == 1:1
        @test colsupport(B,3) == rowsupport(B,3) == 3:3
        @test colsupport(B,5) == rowsupport(B,5) == 5:5
        B = BroadcastArray(*, Diagonal(randn(5)), 2)
        @test B == broadcast(*,B.args...)
        @test colsupport(B,1) == rowsupport(B,1) == 1:1
        @test colsupport(B,3) == rowsupport(B,3) == 3:3
        @test colsupport(B,5) == rowsupport(B,5) == 5:5
        B = BroadcastArray(*, Diagonal(randn(5)), randn(5))
        @test B == broadcast(*,B.args...)
        @test colsupport(B,1) == rowsupport(B,1) == 1:1
        @test colsupport(B,3) == rowsupport(B,3) == 3:3
        @test colsupport(B,5) == rowsupport(B,5) == 5:5

        B = BroadcastArray(+, Diagonal(randn(5)), 2)
        @test colsupport(B,1) == rowsupport(B,1) == 1:5
        @test colsupport(B,3) == rowsupport(B,3) == 1:5
        @test colsupport(B,5) == rowsupport(B,5) == 1:5
    end

    @testset "vector*matrix broadcasting #27" begin
        H = [1., 0.]
        @test Mul(H, H') .+ 1 == H*H' .+ 1
        B =  randn(2,2)
        @test Mul(H, H') .+ B == H*H' .+ B
    end

    @testset "BroadcastArray +" begin
        a = BroadcastArray(+, randn(400), randn(400))
        b = similar(a)
        copyto!(b, a)
        if VERSION ≥ v"1.1"
            @test @allocated(copyto!(b, a)) == 0
        end
        @test b == a
    end

    @testset "Lazy range" begin
        @test broadcasted(LazyArrayStyle{1}(), +, 1:5) ≡ 1:5
        @test broadcasted(LazyArrayStyle{1}(), +, 1, 1:5) ≡ 2:6
        @test broadcasted(LazyArrayStyle{1}(), +, 1:5, 1) ≡ 2:6

        @test broadcasted(LazyArrayStyle{1}(), +, Fill(2,5)) ≡ Fill(2,5)
        @test broadcasted(LazyArrayStyle{1}(), +, 1, Fill(2,5)) ≡ Fill(3,5)
        @test broadcasted(LazyArrayStyle{1}(), +, Fill(2,5), 1) ≡ Fill(3,5)
        @test broadcasted(LazyArrayStyle{1}(), +, Ref(1), Fill(2,5)) ≡ Fill(3,5)
        @test broadcasted(LazyArrayStyle{1}(), +, Fill(2,5), Ref(1)) ≡ Fill(3,5)
        @test broadcasted(LazyArrayStyle{1}(), +, 1, Fill(2,5)) ≡ Fill(3,5)
        @test broadcasted(LazyArrayStyle{1}(), +, Fill(2,5), Fill(3,5)) ≡ Fill(5,5)

        @test broadcasted(LazyArrayStyle{1}(), *, Zeros(5), Zeros(5)) ≡ Zeros(5)
        b = BroadcastArray(exp, randn(5))
        @test b .* Zeros(5) ≡ Zeros(5)
        @test Zeros(5) .* b ≡ Zeros(5)
    end

    @testset "Sub-broadcast" begin
        A = BroadcastArray(exp,randn(5,5))
        V = view(A, 1:2,2:3)
        @test MemoryLayout(typeof(V)) isa BroadcastLayout{typeof(exp)}
        @test BroadcastArray(V) == A[1:2,2:3] == Array(A)[1:2,2:3]

        B = BroadcastArray(-, randn(5,5), randn(5))
        V = view(B, 1:2,2:3)
        @test MemoryLayout(typeof(V)) isa BroadcastLayout{typeof(-)}
        @test BroadcastArray(V) == B[1:2,2:3] == Array(B)[1:2,2:3]
    end

    @testset "AdjTrans" begin
        A = BroadcastArray(exp,randn(5,5))
        @test MemoryLayout(typeof(transpose(A))) isa BroadcastLayout{typeof(exp)}
        @test MemoryLayout(typeof(A')) isa BroadcastLayout{typeof(exp)}
        @test BroadcastArray(A') == BroadcastArray(transpose(A)) == A' == Array(A)'

        B = BroadcastArray(-, randn(5,5), randn(5))
        @test MemoryLayout(typeof(transpose(B))) isa BroadcastLayout{typeof(-)}
        @test MemoryLayout(typeof(B')) isa BroadcastLayout{typeof(-)}
        @test BroadcastArray(B') == BroadcastArray(transpose(B)) == B' == Array(B)'

        Vc = view(B', 1:2,1:3)
        Vt = view(transpose(B), 1:2,1:3)
        @test MemoryLayout(typeof(Vc)) isa BroadcastLayout{typeof(-)}
        @test MemoryLayout(typeof(Vt)) isa BroadcastLayout{typeof(-)}
        @test arguments(Vc) == (B.args[1][1:3,1:2]', permutedims(B.args[2][1:3]))
        @test arguments(Vt) == (transpose(B.args[1][1:3,1:2]), permutedims(B.args[2][1:3]))
        @test BroadcastArray(Vc) == BroadcastArray(Vt) == Vc == Vt == (Array(B)')[1:2,1:3]

        Vc = view(B,1:3,1:2)'
        Vt = transpose(view(B,1:3,1:2))
        @test MemoryLayout(typeof(Vc)) isa BroadcastLayout{typeof(-)}
        @test MemoryLayout(typeof(Vt)) isa BroadcastLayout{typeof(-)}
        @test arguments(Vc) == (B.args[1][1:3,1:2]', permutedims(B.args[2][1:3]))
        @test arguments(Vt) == (transpose(B.args[1][1:3,1:2]), permutedims(B.args[2][1:3]))
        @test BroadcastArray(Vc) == BroadcastArray(Vt) == Vc == (Array(B)')[1:2,1:3]
    end
end


text:
#Test different backends of EasyPlotMPL/PyPlot
#-------------------------------------------------------------------------------
module CMDimData_SampleUsage

using CMDimData
using CMDimData.EasyPlot
CMDimData.@includepkg EasyPlotMPL

#==Constants
===============================================================================#
backendtestlist = [:tk, :gtk3, :gtk, :qt, :wx]
demolist = EasyPlot.demofilelist()

#==Helper functions
===============================================================================#
printsep(label, sep="-") = println("\n", label, "\n", repeat(sep, 80))
printheader(label) = printsep(label, "=")

#==Display sample EasyPlot plots on different PyPlot backends
===============================================================================#
printheader("PyPlot backend test")
plot = evalfile(demolist[1])
for backend in backendtestlist
	printsep("Backend $backend...")
	plot.title = "Backend: $backend"
	try
		disp = EasyPlot.GUIDisplay(:PyPlot, backend=backend)
		display(disp, plot)
	catch e
		@warn e.msg
	end
end

end #module
:SampleCode_Executed


text:
@inline Base.:^(v::AbstractSIMD{W,T}, i::IntegerTypesHW) where {W,T} = Base.power_by_squaring(v, i)
@inline Base.:^(v::AbstractSIMD{W,T}, i::IntegerTypesHW) where {W,T<:Union{Float32,Float64}} = Base.power_by_squaring(v, i)
@inline relu(x) = (y = zero(x); ifelse(x > y, x, y))

@inline Base.fld(x::AbstractSIMD, y::AbstractSIMD) = div(promote_div(x,y)..., RoundDown)
@inline Base.fld(x::AbstractSIMD, y::Real) = div(promote_div(x,y)..., RoundDown)
@inline Base.fld(x::Real, y::AbstractSIMD) = div(promote_div(x,y)..., RoundDown)

@inline function Base.div(x::AbstractSIMD{W,T}, y::AbstractSIMD{W,T}, ::RoundingMode{:Down}) where {W,T<:Integer}
    d = div(x, y)
    d - (signbit(x ⊻ y) & (d * y != x))
end

@inline Base.mod(x::AbstractSIMD{W,T}, y::AbstractSIMD{W,T}) where {W,T<:Integer} =
    ifelse(y == -1, zero(x), x - fld(x, y) * y)

@inline Base.mod(x::AbstractSIMD{W,T}, y::AbstractSIMD{W,T}) where {W,T<:Unsigned} =
    rem(x, y)

@inline function Base.mod(x::AbstractSIMD{W,T1}, y::T2) where {W,T1<:SignedHW,T2<:UnsignedHW}
    _x, _y = promote_div(x, y)
    unsigned(mod(_x, _y))
end
@inline function Base.mod(x::AbstractSIMD{W,T1}, y::T2) where {W,T1<:UnsignedHW,T2<:SignedHW}
    _x, _y = promote_div(x, y)
    signed(mod(_x, _y))
end
@inline function Base.mod(x::AbstractSIMD{W,T1}, y::AbstractSIMD{W,T2}) where {W,T1<:SignedHW,T2<:UnsignedHW}
    _x, _y = promote_div(x, y)
    unsigned(mod(_x, _y))
end
@inline function Base.mod(x::AbstractSIMD{W,T1}, y::AbstractSIMD{W,T2}) where {W,T1<:UnsignedHW,T2<:SignedHW}
    _x, _y = promote_div(x, y)
    signed(mod(_x, _y))
end

@inline Base.mod(i::AbstractSIMD{<:Any,<:Integer}, r::AbstractUnitRange{<:Integer}) =
    mod(i-first(r), length(r)) + first(r)

@inline Base.mod(x::AbstractSIMD, y::NativeTypes) = mod(promote_div(x,y)...)
@inline Base.mod(x::NativeTypes, y::AbstractSIMD) = mod(promote_div(x,y)...)

# avoid ambiguity with clamp(::Missing, lo, hi) in Base.Math at math.jl:1258
# but who knows what would happen if you called it
for (X,L,H) in Iterators.product(fill([:Any, :Missing, :AbstractSIMD], 3)...)
    any(==(:AbstractSIMD), (X,L,H)) || continue
    @eval @inline function Base.clamp(x::$X, lo::$L, hi::$H)
        x_, lo_, hi_ = promote(x, lo, hi)
        ifelse(x_ > hi_, hi_, ifelse(x_ < lo_, lo_, x_))
    end
end

@inline Base.clamp(x::AbstractSIMD{<:Any,<:Integer}, r::AbstractUnitRange{<:Integer}) =
    clamp(x, first(r), last(r))


text:
function HermiteInterpolator()
    return HermiteInterpolator(())
end

function add_sample_point(obj::HermiteInterpolator, arg0::jdouble, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "addSamplePoint", void, (jdouble, Vector{Vector{jdouble}}), arg0, arg1)
end

function derivatives(obj::HermiteInterpolator, arg0::jdouble, arg1::jint)
    return jcall(obj, "derivatives", Vector{Vector{jdouble}}, (jdouble, jint), arg0, arg1)
end

function get_polynomials(obj::HermiteInterpolator)
    return jcall(obj, "getPolynomials", Vector{PolynomialFunction}, ())
end

function value(obj::HermiteInterpolator, arg0::Derivative)
    return jcall(obj, "value", Vector{Derivative}, (Derivative,), arg0)
end

function value(obj::HermiteInterpolator, arg0::DerivativeStructure)
    return jcall(obj, "value", Vector{DerivativeStructure}, (DerivativeStructure,), arg0)
end

function value(obj::HermiteInterpolator, arg0::jdouble)
    return jcall(obj, "value", Vector{jdouble}, (jdouble,), arg0)
end


text:
include("../fileFind/filesfind.jl")
include("../RawFileofflineProc/ReadFiberData.jl")
# using Base.write
include("baseTransferLib.jl")
###########################################################

# endtime = ""
function transFiles(rootdir,startTIme, endtime )
    println()
    files = findAllFiles(rootdir,".conf",startTIme,endtime)
    for t in files
        println( t )
    end
    confgfilename = files[1]
    println( "Selected config file: ",confgfilename )
    println()

    files = findAllFiles(rootdir,".energy",startTIme,endtime)
    transferfiles( files,confgfilename,".RAW4" )

    files = findAllFiles(rootdir,".LC",startTIme,endtime)
    transferfiles( files,confgfilename,".RAW3" )

    println()
    files = findAllFiles(rootdir,".RAW",startTIme,endtime)
    transferfiles( files,confgfilename,".RAW1" )
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
module Yields

import BSplineKit
import ForwardDiff
using LinearAlgebra
using UnicodePlots
using Roots

# don't export type, as the API of Yields.Zero is nicer and
# less polluting than Zero and less/equally verbose as ZeroYieldCurve or ZeroCurve
export rate, discount, accumulation, forward,
LinearSpline, QuadraticSpline, Periodic, Continuous

include("AbstractYield.jl")
include("Rate.jl")
const DEFAULT_COMPOUNDING = Yields.Continuous()

include("utils.jl")
include("bootstrap.jl")
include("SmithWilson.jl")
include("generics.jl")
include("RateCombination.jl")

end


text:
# TODO: temporary until the start of JuliaCon
# JuliaCon.debugmode(true)

function commander(c::Client, m::Message, ::Val{:julia_con})
    startswith(m.content, COMMAND_PREFIX * "con") || return

    regex = Regex(COMMAND_PREFIX * raw"con( help| juliacon2021| now| today)? *(.*)$")

    matches = match(regex, m.content)

    if matches === nothing || matches.captures[1] ∈ (" help", nothing)
        help_commander(c, m, :julia_con)
    elseif matches.captures[1] == " juliacon2021"
        @info "juliacon2021 was required" m.content m.author.username m.author.discriminator
        jc_juliacon2021(c, m)
    elseif matches.captures[1] == " now"
        @info "now was required" m.content m.author.username m.author.discriminator
        jc_now(c, m)
    elseif matches.captures[1] == " today"
        @info "today was required" m.content m.author.username m.author.discriminator
        jc_today(c, m)
    else
        reply(c, m, "Sorry, are you looking for information on JuliaCon 2021? Please check out `con help`")
    end
    return nothing
end

function help_commander(c::Client, m::Message, ::Val{:julia_con})
    reply(c, m, """
        How to look for information about JuliaCon 2021 with the `con` command:
        ```
        con help
        con juliacon2021
        con now
        con today
        ```
        `help` prints this message
        `juliacon2021` greets the users and provide a link to the official website
        `now` details the talks and events occuring right now on the different tracks
        `today` lists the talks and event of the day
        """)
end

function jc_juliacon2021(c::Client, m::Message)
    reply(c, m, "Welcome to JuliaCon 2021! Find more information at https://juliacon.org/2021/.")
    return nothing
end

function jc_now(c::Client, m::Message)
    if JuliaCon.get_running_talks(; now=FAKENOW) === nothing
        @info "There is no JuliaCon program today!"
    else
        reply(c, m, JuliaCon.now(now=FAKENOW, output=:text))
    end
    return nothing
end

function split_pretty_table(str)
    lines = split(str, "\n")
    acc = "```\n"
    strings = Vector{String}()
    for line in lines
        if sum(length, [acc, line]) ≤ 1996
            acc *= "\n" * line
        else
            push!(strings, acc * "\n```")
            acc = "```\n" * line
        end
    end
    push!(strings, acc * "\n```")
    return strings
end

function jc_today(c::Client, m::Message)
    if JuliaCon.get_today(; now=FAKENOW) === nothing
        @info "There is no JuliaCon program today!"
    else
        strings = Vector{String}()
        aux = JuliaCon.today(now = FAKENOW, output = :text)
        tables, legend = aux[1:end-1], aux[end]
        for t in tables, str in split_pretty_table(t)
            push!(strings, str)
        end
        push!(strings, legend)
        foreach(s -> reply(c, m, s), strings)
    end
    return nothing
end


text:
using OrdinaryDiffEq
using DiffEqCallbacks
using RecipesBase
using SimulationLogs
using Suppressor
using Test

function lorenz!(du, u, p, t)
    @log a = u[2]-u[1]
    @log b u[3] + a

    @log c du[1] = p[1]*a
    @log p1 p[1]
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]

    accum = 0
    for i in 1:5
        accum += u[1]
        @log d accum
        if t>50
            @log e u[2]
        end
    end

    return nothing
end

p = [10.0, 28.0, 8/3]
u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 100.0)

prob = ODEProblem(lorenz!, u0, tspan, p)

@testset "Lorenz example" begin
    sol = solve(prob, Tsit5())

    out = get_log(sol)

    a = [u[2]-u[1] for u in sol.u]
    b = [u[3]+(u[2]-u[1]) for u in sol.u]
    c = [p[1]*(u[2]-u[1]) for u in sol.u]
    d = [n*u[1] for u in sol.u, n in 1:5]

    @testset "Get" begin
        @test out.a == a
        @test out.b == b
        @test out.c == c
        @test out.d == d
        @test all(out.p1 .== p[1])

        @test out.a == out[:a]
    end

    @testset "Variable types" begin
        @test typeof(out.a) === typeof(a)
        @test typeof(out.b) === typeof(b)
        @test typeof(out.c) === typeof(c)
        @test typeof(out.d) === typeof(d)
    end

    @testset "Plot recipes" begin
        using SimulationLogs: Scope
        using RecipesBase: apply_recipe

        AnyDict = Dict{Symbol, Any}
        t = range(tspan..., length=2000)
        out = get_log(sol, t)

        rec = apply_recipe(AnyDict(), Scope((sol, :a)))
        @test only(rec).plotattributes == AnyDict(:xguide=>"t", :label=>"a", :seriestype=>:path)
        @test only(rec).args == (t, out.a)

        rec = apply_recipe(AnyDict(), Scope((sol, [:b, :c])))
        @test rec[1].plotattributes == AnyDict(:xguide=>"t", :label=>"b", :seriestype=>:path)
        @test rec[1].args == (t, out.b)
        @test rec[2].plotattributes == AnyDict(:xguide=>"t", :label=>"c", :seriestype=>:path)
        @test rec[2].args == (t, out.c)

        rec = apply_recipe(AnyDict(), Scope((sol, (:a, :c))))
        @test only(rec).plotattributes == AnyDict(:xguide=>"t", :label=>"(a, c)", :seriestype=>:path)
        @test only(rec).args == (out.a, out.c)

        rec = apply_recipe(AnyDict(), Scope((sol, :d)))
        @test only(rec).plotattributes == AnyDict(:xguide=>"t", :label=>["d[:,$i]" for i in (1:5)'], :seriestype=>:path)
        @test only(rec).args == (t, out.d)
    end

    @testset "Printing" begin
        @test @capture_out(print(out)) ==
        """
        SimulationLog with signals:
          a :: $(typeof(out.a))
          b :: $(typeof(out.b))
          d :: $(typeof(out.d))
          c :: $(typeof(out.c))
          p1 :: $(typeof(out.p1))
        """
    end

    @testset "Miscellaneous" begin
        @test Set(propertynames(out)) == Set((:a, :b, :c, :d, :p1))

        @test is_active(out) === false

        @test_logs () get_log(sol)
    end
end

@testset "Callbacks" begin
    d_callback = PresetTimeCallback(Int(tspan[1]):Int(tspan[2]), integrator->(integrator.p[1] += 0.01))
    callback_set = CallbackSet(d_callback)
    prob.p[1] = 10.0
    sol = logged_solve(prob, Tsit5(); callback=callback_set)
    prob.p[1] = 10.0
    out = get_log(sol)

    @test out.c[1] == p[1]*(u0[2]-u0[1])
    @test !all(out.c .≈ [p[1]*(u[2]-u[1]) for u in sol.u])
    @test !all(==(p[1]), out.p1)
    @test out.p1[end] ≈ 11
end

@testset "Logged solution" begin
    include("lotka_test.jl")
end

@testset "Global log" begin
    using SimulationLogs: activate!, deactivate!, reset!

    glog = SimulationLogs.GLOBAL_LOG

    @test propertynames(glog) === ()

    @test collect(keys(glog)) == Symbol[]

    @test is_active(glog) === false
    activate!()
    @test is_active(glog) === true
    deactivate!()

    @test isempty(value_dict(glog))
    glog[:a] = zeros(5)
    @test glog.a == zeros(5)
    @test glog[:a] == zeros(5)
    reset!()
    @test_throws KeyError glog.a
    @test_throws KeyError glog[:a]
    @test isempty(value_dict(glog))
end

max_stars_repo_name: jonniedie/SimulationLogs

max_issues_repo_name: jonniedie/SimulationLogs

max_forks_repo_name: jonniedie/SimulationLogs


text:
# -------------------------------------------------------------------------------------
# parse_arg
"""
This is part of the public interface.
Implement a new method for parsing a new type.
`parse_arg` must return one of the following types
[`Symbol`, `Int`, `Float64`, `String`, `Bool`, `NamedTuple`, `Dict`]
or other type that already has an implemented `parse_arg` method.
For the case in `NamedTuple` and `Dict`, the keys and values must be
one of the first other types.

# Example:
```julia
struct Foo
    x::Int
    y::String
end
parse_arg(f::Foo) = (;f.x, f.y)

foo = Foo(1, "A")
fname = dfname("my_foo", foo, ".ext")
# "my_foo [x=1 y=A].ext"
```
"""
parse_arg
parse_arg(s::Symbol) = s
parse_arg(n::Integer) = n
parse_arg(f::AbstractFloat) = f
parse_arg(s::AbstractString) = string(s)
parse_arg(b::Bool) = b
parse_arg(p::Pair) = p
parse_arg(pt::NamedTuple) = pt
parse_arg(pt::Iterators.Pairs) = parse_arg(Dict(pt...))
parse_arg(pt::Dict) = pt
parse_arg(v::Any) = error("parse_arg(v::", typeof(v), ") not implemented. Type `?parse_arg` for help.")

# -------------------------------------------------------------------------------------
# _argstr
_argstr(s::Symbol) = _check_str(string(s))
_argstr(n::Integer) = _check_str(string(n))
_argstr(f::AbstractFloat) = _check_str(@sprintf("%0.2e", f))
_argstr(s::AbstractString) = _check_str(string(s))
_argstr(b::Bool) = _check_str(string(b))

## ----------------------------------------------------------------------------
# Will extract the first Vector{<:AbstractString} args and makes
# them a path and returns the rest
function _extract_dir(args...)
    isempty(args) && return ("", args)
    path = ""
    for (i, arg) in enumerate(args)
        !(arg isa Vector{<:AbstractString}) && return (path, args[i:end])
        for dir in arg
            path = joinpath(path, dir)
        end
    end
    return (path, tuple())
end

## ----------------------------------------------------------------------------
# assumes _extract_dir was already called
# extract all first _isvalT ignoring the last arg
function _extract_head(args...)
    head_args = []
    isempty(args) && return (head_args, args)
    lasti = lastindex(args)
    for (i, arg) in enumerate(args)
        parg = parse_arg(arg)
        if ((i != lasti) && _isvalT(parg)); push!(head_args, parg)
            else; return (head_args, args[i:end])
        end
    end
    return (head_args, tuple(args[lasti]))
end

# -------------------------------------------------------------------------------------
# assumes _extract_head was already called
# extract all first _ispairT
function _extract_params(args...)
    params_args = []
    isempty(args) && return (params_args, args)
    for (i, arg) in enumerate(args)
        parg = parse_arg(arg)
        if _ispairT(parg); push!(params_args, parg)
            else; return (params_args, args[i:end])
        end
    end
    return (params_args, tuple())
end

# -------------------------------------------------------------------------------------
function _extract_ext(args...)
    len = length(args)
    for arg in args
        parg = parse_arg(arg)
        _isT = _isvalT(parg)
        # valid ext arg
        (_isT && len == 1) && return parg
        # several args
        _isT && error(
            "After the first key:value argument no single value is allowed (except the extension at the end). ",
            "Type protocole_desc() for info!"
        )
        # here such arg should also fails _ispairT
        error("parse_arg(a::$(typeof(arg))) returns an invalid type. Use `?parse_arg` for help.")
    end
    return "" # default ext
end

# -------------------------------------------------------------------------------------
# dfname
function _dfname(args...)
    _check__SEPS()

    # --------------------------------------------------------
    # extract args
    head_args, args = _extract_head(args...)
    params_args, args = _extract_params(args...)

    # --------------------------------------------------------
    # deal with ext
    extarg = _extract_ext(args...)
    sextarg = _argstr(extarg)
    if startswith(sextarg, _SEPS[:EXT_SEP]) || isempty(sextarg)
        ext = sextarg
    elseif length(params_args) > 0
        ext = string(_SEPS[:EXT_SEP], sextarg)
    else
        # is not an extension
        ext = _argstr("")
        push!(head_args, sextarg)
    end

    # --------------------------------------------------------
    # join params to check for collisions
    params_dict = Dict()
    for param_arg in params_args
        param_arg = (param_arg isa Pair) ? [param_arg] : pairs(param_arg)
        for (k, v) in param_arg
            ck = _check_keyT(k)
            cv = _check_valT(v)
            strk = _argstr(ck)
            strv = _argstr(cv)

            haskey(params_dict, strk) &&
                error("You have passed two keys that lead to the same string '", strk, "'. Collisions are not allowed")
            params_dict[strk] = strv
        end
    end

    # --------------------------------------------------------
    # build name
    head_strs = [_argstr(arg) for arg in head_args]
    params_scol = sort!(collect(params_dict); by = first)
    params_strs = [string(k, _SEPS[:PAIR_SEP], v) for (k, v) in params_scol]
    if !isempty(params_strs)
        param_str = string(
            _SEPS[:PARAMS_LSEP],
            join(params_strs, _SEPS[:ELEMT_SEP]),
            _SEPS[:PARAMS_RSEP]
        )
        push!(head_strs, param_str)
    end
    filter!(!isempty, head_strs)
    fname = string(join(head_strs, _SEPS[:ELEMT_SEP]), ext)

    # --------------------------------------------------------
    # happyness
    return fname
end

_dfname() = ""
_dfname(fname::String) = isvalid_dfname(fname) ? fname : _dfname("", fname)

function dfname(args...)
    dir, args = _extract_dir(args...)
    fname = _dfname(args...)
    isempty(fname) ? dir : joinpath(dir, fname)
end


text:
@testset "AlgorithmConfig" begin
    let config = Cerberus.AlgorithmConfig()
        @test config.lp_solver_factory == Cerberus.DEFAULT_LP_SOLVER_FACTORY
        @test config.silent == Cerberus.DEFAULT_SILENT
        @test config.branching_rule == Cerberus.DEFAULT_BRANCHING_RULE
        @test config.time_limit_sec == Cerberus.DEFAULT_TIME_LIMIT_SEC
        @test config.node_limit == Cerberus.DEFAULT_NODE_LIMIT
        @test config.gap_tol == Cerberus.DEFAULT_GAP_TOL
        @test config.int_tol == Cerberus.DEFAULT_INTEGRALITY_TOL
        @test config.warm_start_strategy == Cerberus.DEFAULT_WARM_START_STRATEGY
        @test config.model_reuse_strategy ==
              Cerberus.DEFAULT_MODEL_REUSE_STRATEGY
        @test config.formulation_tightening_strategy ==
              Cerberus.DEFAULT_FORMULATION_TIGHTENING_STRATEGY
    end

    let lp = identity,
        sl = true,
        br = Cerberus.PseudocostBranching(),
        tl = 12.3,
        nl = 10,
        gt = 10.0,
        it = 1e-6,
        ws = Cerberus.NO_WARM_STARTS,
        mr = Cerberus.NO_MODEL_REUSE,
        ts = Cerberus.STATIC_FORMULATION,

        config = Cerberus.AlgorithmConfig(
            lp_solver_factory = lp,
            silent = sl,
            branching_rule = br,
            time_limit_sec = tl,
            node_limit = nl,
            gap_tol = gt,
            int_tol = it,
            warm_start_strategy = ws,
            model_reuse_strategy = mr,
            formulation_tightening_strategy = ts,
        )

        @test config.lp_solver_factory == lp
        @test config.silent == sl
        @test config.branching_rule == br
        @test config.time_limit_sec == tl
        @test config.node_limit == nl
        @test config.gap_tol == gt
        @test config.int_tol == it
        @test config.warm_start_strategy == ws
        @test config.model_reuse_strategy == mr
        @test config.formulation_tightening_strategy == ts
    end
end


text:
using Test

mutable struct MyType
   arr::Array{Float64,1}
   var::Int
   str::String
   MyType(x...) = new(x[1],x[2],x[3])
   MyType() = new()
end

function f(M::MyType)
   mm = MyType()
   mm.arr = M.arr
   return mm
end

@testset "Ensure memory return is working properly..." begin

m = MyType(rand(3),2,"hello")

mmm = f(m)

mmm.arr[1] = 1.0

@test m.arr[1] == 1.0

end


text:
# # DIRS
# const MODEL_RAW_DIR = joinpath(RAW_DATA_DIR, PROJ_IDER)
# const MODEL_PROCESSED_DATA_DIR = joinpath(PROCESSED_DATA_DIR, basename(MODEL_RAW_DIR))
# const MODEL_FIGURES_DIR = joinpath(FIGURES_DATA_DIR, basename(MODEL_RAW_DIR))
# const MODEL_CACHE_DIR = joinpath(MODEL_PROCESSED_DATA_DIR, "cache")

# function _create_dirs()
#     for dir in [MODEL_PROCESSED_DATA_DIR,
#                 MODEL_FIGURES_DIR,
#                 MODEL_CACHE_DIR
#             ]
#         isdir(dir) && continue
#         try; mkpath(dir); catch end
#     end
# end

# # FILES
# const MODEL_RAW_MAT_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "iJR904.mat")
# const BASE_MODELS_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "base_models.bson")
# const EXCH_MET_MAP_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "exch_met_map.bson")
# const MAXENT_B0SEEDS_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "maxent_b0seed_bundles.bson")
# const MAXENT_VARIANTS_INDEX_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "maxent_variants_index.bson")
# const CORR_DAT_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "corr_dat_file.bson")
# const LP_DAT_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "lp_dat.bson")


text:
export optimize

"""
$(TYPEDSIGNATURES)

Optimizes the `HermiteMapComponent` `C` based on the ensemble matrix `X`. Several kind of optimization are implemented
depending on the desired robustness of the map

Several options for `optimkind` are available:
* `kfold` uses a k-fold cross validation procedure (the more robust choice)
* `split` splits the set of samples into a training and a testing
*  An `Int64` to determine the maximum number of features
* `nothing` to simply optimize the existing coefficients in the basis expansion.

The following optionial settings can be tuned:
* `maxterms::Int64 = 1000`: a maximum number of terms for each map component
* `withconstant::Bool = false`: the option to remove the constant feature in the greedy optimization routine, if the zero feature is the constant function for the basis of interest.
* `withqr::Bool = false`: improve the conditioning of the optimization problem with a QR factorization of the feature basis (recommended option)
* `maxpatience::Int64 = 10^5`: for `optimkind = split`, the maximum number of extra terms that can be added without decreasing the validation error before the greedy optimmization get stopped.
* `verbose::Bool = false`: prints details of the optimization procedure, current component optimize, training and validation errors, number of features added.
* `hessprecond::Bool=true`: use a preconditioner based on the Gauss-Newton of the Hessian of the loss function to accelerate the convergence.
* `start::Int64=1`: the first component of the map to optimize
* `P::Parallel = serial`: option to use multi-threading to optimize in parallel the different map components
* `ATMcriterion::String="gradient"`: sensitivty criterion used to select the feature to add to the expansion. The default uses the derivative of the cost function with respect to the coefficient
   of the features in the reduced margin of the current set of features.
"""
function optimize(C::HermiteMapComponent, X, optimkind::Union{Nothing, Int64, String};
                  maxterms::Int64 = 100,
                  withconstant::Bool = false, withqr::Bool = false,
                  maxpatience::Int64 = 10^5, verbose::Bool = false,
                  hessprecond = true, P::Parallel = serial, ATMcriterion::String = "gradient")

    m = C.m
    Nx = C.Nx

    if verbose == true
        println("Optimizing component "*string(Nx)*":")
    end

    if typeof(optimkind) <: Nothing
        S = Storage(C.I.f, X)

        # Optimize coefficients
        if withqr == false
            coeff0 = getcoeff(C)

            if hessprecond == true
                precond = zeros(ncoeff(C), ncoeff(C))
                precond!(precond, coeff0, S, C, X)
                precond_chol = cholesky(Symmetric(precond); check = false)

                if issuccess(precond_chol) == true
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10, P = Preconditioner(Symmetric(precond), precond_chol)))
                elseif cond(Diagonal(precond)) < 10^6  #try the diagonal preconditioner
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10), P = Diagonal(precond))
                else # don't use any preconditioner
                    res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                          Optim.LBFGS(; m = 10))
                end
            else
                res = Optim.optimize(Optim.only_fg!(negative_log_likelihood(S, C, X)), coeff0,
                      Optim.LBFGS(; m = 10))
            end

            setcoeff!(C, Optim.minimizer(res))
            error = res.minimum
        else
            F = QRscaling(S)
            coeff0 = getcoeff(C)
            mul!(coeff0, F.U, coeff0)

            # mul!(S.ψoffψd0, S.ψoffψd0, F.Uinv)
            # mul!(S.ψoffdψxd, S.ψoffdψxd, F.Uinv)
            if hessprecond == true
                qrprecond = zeros(ncoeff(C), ncoeff(C))
                qrprecond!(qrprecond, coeff0, F, S, C, X)
                qrprecond_chol = cholesky(Symmetric(qrprecond); check = false)

                if issuccess(qrprecond_chol) == true
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10, P = Preconditioner(Symmetric(qrprecond), qrprecond_chol)))
                elseif cond(Diagonal(qrprecond_chol)) < 10^6
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10, P = Diagonal(qrprecond)))
                else
                    res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                         Optim.LBFGS(; m = 10))
                end
            else
                res = Optim.optimize(Optim.only_fg!(qrnegative_log_likelihood(F, S, C, X)), coeff0,
                                     Optim.LBFGS(; m = 10))
            end

            if Optim.converged(res)
                mul!(view(C.I.f.coeff,:), F.Uinv, Optim.minimizer(res))
            else
                error("Optimization hasn't converged")
            end

            error = res.minimum
        end

    elseif typeof(optimkind) <: Int64
        C, error =  greedyfit(m, Nx, X, optimkind;
                              α = C.α, withconstant = withconstant,
                              withqr = withqr, maxpatience = maxpatience,
                              verbose = verbose, hessprecond = hessprecond, b = getbasis(C), ATMcriterion = ATMcriterion)

    elseif optimkind ∈ ("kfold", "kfolds", "Kfold", "Kfolds")
        # Define cross-validation splits of data
        n_folds = 5
        folds = kfolds(1:size(X,2), k = n_folds)

        # Run greedy approximation
        max_iter = min(m-1, ceil(Int64, sqrt(size(X,2))), maxterms)

        valid_error = zeros(max_iter, n_folds)
        if typeof(P) <: Serial
            @inbounds for i=1:n_folds
                idx_train, idx_valid = folds[i]

                if verbose == true
                    println("Fold "*string(i)*" / "*string(n_folds)*":")
                end
                C, error = greedyfit(m, Nx, X[:,idx_train], X[:,idx_valid], max_iter;
                                     α = C.α, withconstant = withconstant, withqr = withqr, verbose  = verbose,
                                     hessprecond = hessprecond, b = getbasis(C), ATMcriterion = ATMcriterion)

                # error[2] contains the history of the validation error
                valid_error[:,i] .= deepcopy(error[2])
            end
        elseif typeof(P) <: Thread
            @inbounds  Threads.@threads for i=1:n_folds
                idx_train, idx_valid = folds[i]

                if verbose == true
                    println("Fold "*string(i)*":")
                end

                C, error = greedyfit(m, Nx, X[:,idx_train], X[:,idx_valid], max_iter;
                                     α = C.α, withconstant = withconstant, withqr = withqr, verbose  = verbose,
                                     hessprecond = hessprecond, b = getbasis(C), ATMcriterion = ATMcriterion)

                # error[2] contains the history of the validation error
                valid_error[:,i] .= deepcopy(error[2])
            end
        end

        # Find optimal numbers of terms
        mean_valid_error = mean(valid_error, dims  = 2)[:,1]

        _, opt_nterms = findmin(mean_valid_error)

        # Run greedy fit up to opt_nterms on all the data
        C, error = greedyfit(m, Nx, X, opt_nterms;
                             α = C.α, withqr = withqr, verbose  = verbose,
                             hessprecond = hessprecond, b = getbasis(C), ATMcriterion = ATMcriterion)

    elseif optimkind ∈ ("split", "Split")
        # 20% of the data is used for the cross-validation
        valid_train_split = 0.2
        nvalid = ceil(Int64, floor(valid_train_split*size(X,2)))
        X_train = X[:,nvalid+1:end]
        X_valid = X[:,1:nvalid]

        # Run greedy approximation
        max_iter =  min(m-1, maxterms, ceil(Int64, sqrt(size(X,2))))

        Cvalid, error = greedyfit(m, Nx, X_train, X_valid, max_iter;
                             α = C.α, withconstant = withconstant, withqr = withqr,
                             maxpatience = maxpatience, verbose  = verbose,
                             hessprecond = hessprecond, b = getbasis(C), ATMcriterion = ATMcriterion)

        # Find optimal numbers of terms
        train_error, valid_error = error
        _, opt_nterms = findmin(valid_error)

        # Since we add the features ina greedy fashion, we can simply pick the number of features that minimize the loss functino,
        # and optimize the coefficients with the entire data set

        # Run greedy fit up to opt_nterms on all the data
        C = HermiteMapComponent(Cvalid.m, Cvalid.Nx, getidx(Cvalid)[1:opt_nterms, :],
                                zeros(opt_nterms); α = Cvalid.α, b = getbasis(Cvalid))

        C, error = optimize(C, X, nothing;
                 maxterms = maxterms, withconstant = withconstant, withqr = withqr,
                 maxpatience = maxpatience, verbose = false, hessprecond = hessprecond,
                 P = P, ATMcriterion = ATMcriterion)
    else
        error("Argument max_terms is not recognized")
    end
    return C, error
end

"""
$(TYPEDSIGNATURES)

Optimizes the `LinHermiteMapComponent` `L` based on the ensemble matrix `X`. Several kind of optimization are implemented
depending on the desired robustness of the map

Several options for `optimkind` are available:
* `kfold` uses a k-fold cross validation procedure (the more robust choice)
* `split` splits the set of samples into a training and a testing
*  An `Int64` to determine the maximum number of features
* `nothing` to simply optimize the existing coefficients in the basis expansion.

The following optionial settings can be tuned:
* `maxterms::Int64 = 1000`: a maximum number of terms for each map component
* `withconstant::Bool = false`: the option to remove the constant feature in the greedy optimization routine, if the zero feature is the constant function for the basis of interest.
* `withqr::Bool = false`: improve the conditioning of the optimization problem with a QR factorization of the feature basis (recommended option)
* `maxpatience::Int64 = 10^5`: for `optimkind = split`, the maximum number of extra terms that can be added without decreasing the validation error before the greedy optimmization get stopped.
* `verbose::Bool = false`: prints details of the optimization procedure, current component optimize, training and validation errors, number of features added.
* `hessprecond::Bool=true`: use a preconditioner based on the Gauss-Newton of the Hessian of the loss function to accelerate the convergence.
* `start::Int64=1`: the first component of the map to optimize
* `P::Parallel = serial`: option to use multi-threading to optimize in parallel the different map components
* `ATMcriterion::String="gradient"`: sensitivty criterion used to select the feature to add to the expansion. The default uses the derivative of the cost function with respect to the coefficient
   of the features in the reduced margin of the current set of features.
"""
function optimize(L::LinHermiteMapComponent, X::Array{Float64,2}, optimkind::Union{Nothing, Int64, String};
                  withconstant::Bool = false, withqr::Bool = false, maxpatience::Int64=20, verbose::Bool = false,
                  hessprecond::Bool = true, ATMcriterion::String="gradient")

    transform!(L.L, X)
    C = L.C
    C_opt, error = optimize(C, X, optimkind; withconstant = withconstant, withqr = withqr, maxpatience = maxpatience,
                            verbose = verbose, hessprecond = hessprecond, ATMcriterion = ATMcriterion)

    itransform!(L.L, X)

    return LinHermiteMapComponent(L.L, C_opt), error
end


text:
using ThreeBodyProblem
using DifferentialEquations
using Plots
gr()

sys = jupiter_europa()

rv₀ = [1.000045941880842, -0.000083429669677, 0.002323145484474, -0.093200789960739, 0.099227116703640, -0.033175567182715]

tspan = (0.,-22.)
prob = ODEProblem(CR3BPdynamics!,rv₀,tspan,sys)
sol = solve(prob,reltol=1e-6)

c = [1-sys.μ, 0, 0]
xyz_hemisphere = deserno_hemisphere(100, c)

plot(sol,vars=(1,2),title="Lunar return trajectory in rotating frame")
plot!(sys)


text:
printstyled("  safe_isfile\n", color=:light_green)
@test safe_isfile("runtests.jl") == true
@test safe_isfile("foo.jl") == false

printstyled("  safe_isdir\n", color=:light_green)
@test safe_isdir("SampleFiles") == true
@test safe_isdir("Roms") == false

printstyled("  ls\n", color=:light_green)

cfile = path*"/SampleFiles/Restricted/03_02_27_20140927.euc.ch"
@test any([occursin("test", i) for i in ls()])

S = [
      "CoreUtils/",
      "CoreUtils/test_ls.jl",
      "CoreUtils/poo",
      "SampleFiles/UW/*W",
      "SampleFiles/UW/02*o",
      "CoreUtils/test_*"
    ]
S_expect =  [
              ["test_calculus.jl", "test_ls.jl", "test_poly.jl", "test_time.jl", "test_typ2code.jl"],
              ["test_ls.jl"],
              String[],
              ["00012502123W", "99011116541W"],
              ["02062915175o", "02062915205o"],
              ["test_calculus.jl", "test_ls.jl", "test_poly.jl", "test_time.jl", "test_typ2code.jl"],
            ]

# Test that ls returns the same files as `ls -1`
for (n,v) in enumerate(S)
  files = String[splitdir(i)[2] for i in ls(v)]
  # if Sys.iswindows() == false
    expected = S_expect[n]
    @test files == expected
  # end
  [@test isfile(f) for f in ls(v)]
end
# Test that ls invokes find_regex under the right circumstances
@test change_sep(ls(S[5])) == change_sep(regex_find("SampleFiles/", r"02.*o$"))

if safe_isfile(cfile)
  T = path .* [
                "/SampleFiles/Restricted/*.cnt",
                "/SampleFiles/*",
                "/SampleFiles/Restricted/2014092709*cnt"
              ]
  T_expect =  [63, 535, 60]
  if safe_isfile(path .* "/SampleFiles/restricted.tar.gz")
    T_expect[2] += 1
  end

  # Test that ls finds the same number of files as bash `ls -1`
  for (n,v) in enumerate(T)
    files = ls(v)
    @test (isempty(files) == false)
    @test (length(files) == T_expect[n])
    [@test isfile(f) for f in files]
  end

  # Test that ls invokes find_regex under the right circumstances
  @test change_sep(ls(T[2])) == change_sep(regex_find("SampleFiles", r".*$"))
  @test change_sep(ls(T[3])) == change_sep(regex_find("SampleFiles", r"Restricted/2014092709.*cnt$"))
else
  printstyled("  extended ls tests skipped. (files not found; is this Appveyor?)\n", color=:green)
end

if Sys.iswindows()
  @test safe_isfile("http://google.com") == false
  @test safe_isdir("http://google.com") == false
end


text:
# estimate the mean by numerical integration over uniform percentiles
estimateMean(d::ContinuousUnivariateDistribution;kwargs...) =
  meanFunOfProb(d;kwargs...,fun=(d,p)->quantile(d,p))

# estimate variance by numerical integration over uniform percentiles
function estimateVariance(d::ContinuousUnivariateDistribution; mean=missing, kwargs...)
    m = ismissing(mean) ? Distributions.mean(d) : mean
    function squaredDiff(d,p)
        t = quantile(d, p) - m
        t*t
    end
    meanFunOfProb(d;kwargs...,fun=squaredDiff)
end

"""
compute mean over a function(uniformly distributed probabilities)

used to estimate moments of logitnorm
"""
function meanFunOfProb(d::ContinuousUnivariateDistribution;relPrec = 1e-4, maxCnt=2^18, fun=(d,p)->quantile.(d,p) )
    δ = 1/32 # start with 31 points (32 intervals between 0 and 1)
    p = δ:δ:(1-δ)
    # for K=1/δ intervals, there are (K-1) points at c_i
    # The first points at δ represents interval (δ/2,3/2δ)
    # The following picture shows points and intervals for K = 4
    #---|---|---|---#
    # |---|---|---| #
    # we need to add points for δ/4 and 1-δ/4 representing the edges
    # but their weight is only half, because they represents half an inverval
    #m = sum(c_i*δ) + el*(δ/2) + er*(δ/2) = (sum(c_i) + er/2 + el/2)*δ
    s = sum(fun.(d,p))   # sum at points c_i
    el = fun(d,δ/4)  #
    er = fun(d,1-δ/4)
    m = (s + el/2 + er/2)*δ
    relErr = 1
    while 1/δ < maxCnt
        mPrev = m
        δ  = δ / 2
        # to double the number of reference points, half the interval
        # for each second point we already computed fun
        # only need to add the new points to the sum of central points
        p = δ:δ*2:(1-δ) # points at the center of current intervals
        s += sum(fun.(d,p))
        el = fun(d,δ/4)
        er = fun(d,1-δ/4)
        m = (s + el/2 + er/2)*δ
        relErr = abs(m - mPrev)/m
        #println("cnt=$(1/δ), m=$m, mPrev=$mPrev, relErr=$relErr")
        #if the estimate did not change much, can return
        relErr <= relPrec && break
    end
    relErr > relPrec && @warn "Returning meanFunOfProb results of low relative precision of $relErr"
    m
end


text:
module Client

# --------------------------------
# Include
# --------------------------------
include("node.jl")
include("app.jl")

# --------------------------------
# Export
# --------------------------------
export Node, register_to_server, fit_service, start

end # module


text:
@testset "reading and writing" begin
    @testset "$fmt - $k" for fmt in (:bson, :julia_serialize), (k, v) in datas
        io = IOBuffer()
        orig = JLSOFile(v; format=fmt)
        write(io, orig)

        seekstart(io)

        result = read(io, JLSOFile)
        @test result.version == orig.version
        @test result.julia == orig.julia
        @test result.image == orig.image
        @test result.manifest == orig.manifest
        @test result.format == orig.format
        @test result.compression == orig.compression
        @test result.objects == orig.objects
        @test result == orig
    end
end

@testset "deserialization" begin
    # Test deserialization works
    @testset "$fmt - $k" for fmt in (:bson, :julia_serialize), (k, v) in datas
        jlso = JLSOFile(v; format=fmt)
        @test jlso["data"] == v
    end

    @testset "unsupported julia version $jlso_version" for jlso_version in (v"1.0", v"2.0")
        jlso = @suppress_out begin
            JLSOFile(
                jlso_version,
                VERSION,
                :julia_serialize,
                :none,
                img,
                pkgs,
                Dict("data" => hw_5),
            )
        end

        # Test failing to deserialize data because of incompatible julia versions
        # will return the raw bytes
        result = if VERSION < v"1.2"
            @test_warn(LOGGER, r"MethodError*", jlso["data"])
        else
            @test_warn(LOGGER, r"TypeError*", jlso["data"])
        end

        @test result == hw_5

        # TODO: Test that BSON works across julia versions using external files?
    end

    @testset "missing module" begin
        # We need to load and use AxisArrays on another process to cause the
        # deserialization error
        pnum = first(addprocs(1))

        try
            # We need to do this separately because there appears to be a race
            # condition on AxisArrays being loaded.
            f = @spawnat pnum begin
                @eval Main using Serialization
                @eval Main using BSON
                @eval Main using AxisArrays
            end

            fetch(f)

            @testset "serialize" begin
                f = @spawnat pnum begin
                    io = IOBuffer()
                    serialize(
                        io,
                        AxisArray(
                            rand(20, 10),
                            Axis{:time}(14010:10:14200),
                            Axis{:id}(1:10)
                        )
                    )
                    return io
                end

                io = fetch(f)
                bytes = take!(io)

                jlso = JLSOFile(
                    v"3.0",
                    VERSION,
                    :julia_serialize,
                    :none,
                    img,
                    project,
                    manifest,
                    Dict(:data => bytes),
                )

                # Test failing to deserailize data because of missing modules will
                # still return the raw bytes
                result = @test_warn(LOGGER, r"KeyError*", jlso[:data])
                @test result == bytes
            end

            @testset "bson" begin
                f = @spawnat pnum begin
                    io = IOBuffer()
                    bson(
                        io,
                        Dict(
                            :data => AxisArray(
                                rand(20, 10),
                                Axis{:time}(14010:10:14200),
                                Axis{:id}(1:10)
                            )
                        )
                    )
                    return io
                end

                io = fetch(f)
                bytes = take!(io)

                jlso = JLSOFile(
                    v"3.0",
                    VERSION,
                    :bson,
                    :none,
                    img,
                    project,
                    manifest,
                    Dict(:data => bytes),
                )

                # Test failing to deserailize data because of missing modules will
                # still return the raw bytes
                result = @test_warn(LOGGER, r"UndefVarError*", jlso["data"])

                @test result == bytes
            end
        finally
            rmprocs(pnum)
        end
    end
end

@testset "saving and loading" begin
    mktempdir() do path
        @testset "$fmt - $k" for fmt in (:bson, :julia_serialize), (k, v) in datas
            JLSO.save("$path/$fmt-$k.jlso", k => v; format=fmt)
            result = JLSO.load("$path/$fmt-$k.jlso")
            @test result[k] == v
        end
    end
end


text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# Example: ELMOD-DE

#using Joulia
include("src/Joulia.jl")

#Run this only once to set up packages
using Pkg
Pkg.activate(pwd())
Pkg.add("DataFrames")
Pkg.add("CSV")
Pkg.add("JuMP")
Pkg.add("ProgressMeter")
Pkg.add(Pkg.PackageSpec(name = "Gurobi", version = v"0.8"))
#Pkg.add("Clp")
#end

using DataFrames
using CSV
using Gurobi
#using Clp

# data load for 2015 sample data
# see http://doi.org/10.5281/zenodo.1044463
pp_df = CSV.read("data_test/power_plants.csv")
avail_con_df = CSV.read("data_test/avail_con.csv")
prices_df = CSV.read("data_test/prices.csv")

storages_df = CSV.read("data_test/storages.csv")

lines_df = CSV.read("data_test/lines.csv")

load_df = CSV.read("data_test/load.csv")
nodes_df = CSV.read("data_test/nodes.csv")
exchange_df = CSV.read("data_test/exchange.csv")

res_df = CSV.read("data_test/res.csv")
avail_pv = CSV.read("data_test/avail_pv.csv")
avail_windon = CSV.read("data_test/avail_windon.csv")
avail_windoff = CSV.read("data_test/avail_windoff.csv")

avail = Dict(:PV => avail_pv,
	:WindOnshore => avail_windon,
	:WindOffshore => avail_windoff,
	:global => avail_con_df)

# generation of Joulia data types
pp = PowerPlants(pp_df, avail=avail_con_df, prices=prices_df)
storages = Storages(storages_df)
lines = Lines(lines_df)
nodes = Nodes(nodes_df, load_df, exchange_df)
res = RenewableEnergySource(res_df, avail)

# generation of the Joulia model
elmod = JouliaModel(pp, res, storages, nodes, lines)

# sclicing the data in weeks with the first full week starting at hour 49
slices = week_slices(49)

# running the Joulia model for week 30 using the Gurobi solver
results = run_model(elmod, slices[30], solver=GurobiSolver())


text:
export compute_mutations!, make_mutdict!, make_mutdict

"""
	compute_mutations(tree::Tree{EvoData})

Compute mutation on each branch of `tree`. Each node of `tree` must have a sequence.
"""
function compute_mutations!(tree::Tree{EvoData})
	for n in values(tree.nodes)
		if !n.isroot
			if 	isempty(n.data.sequence) || isempty(n.anc.data.sequence)
				error("Node $(n.label) or its ancestor does not have a sequence.")
			elseif length(n.data.sequence) != length(n.anc.data.sequence)
				error("Node $(n.label) and its ancestor do not have sequences of the same length")
			else
				if !isempty(n.data.mutations)
					n.data.mutations = Array{Mutation,1}(undef, 0)
				end
				for i in 1:length(n.data.sequence)
					if n.data.sequence[i] != n.anc.data.sequence[i]
						push!(n.data.mutations, Mutation(i, n.anc.data.sequence[i], n.data.sequence[i]))
					end
				end
			end
		end
	end
end

"""
"""
function make_mutdict!(tree::Tree{EvoData}, labellist)
	compute_mutations!(tree)
	mutdict = Dict{Tuple{Int64,Int64,Int64}, Int64}()
	for l in labellist
		for m in tree.lnodes[l].data.mutations
			key = map(f->getfield(m, f), fieldnames(Mutation))
			mutdict[key] = get(mutdict, key, 0) + 1
		end
	end
	return mutdict
end

"""
	make_mutdict!(tree:Tree{EvoData})

Make a dictionary of mutations that appear in `tree`, mapping each mutation to the number of times it appears.
"""
function make_mutdict!(tree::Tree{EvoData}; gaps=false)
	compute_mutations!(tree)
	mutdict = Dict{Tuple{Int64,Int64,Int64}, Int64}()
	mutloc = Dict{Tuple{Int64,Int64,Int64}, Array{String,1}}()
	make_mutdict!(mutdict, mutloc, tree.root, gaps)
	return mutdict, mutloc
end

"""
"""
function make_mutdict!(mutdict, mutloc, node::TreeNode{EvoData}, gaps::Bool)
	for m in node.data.mutations
		key = map(f->getfield(m, f), fieldnames(Mutation))
		if !gaps || key[2]<5 || key[3]<5
			mutdict[key] = get(mutdict, key, 0) + 1
			if length(get(mutloc,key,[])) == 0
				mutloc[key] = [node.label]
			else
				push!(mutloc[key], node.label)
			end
		end
	end
	if !node.isleaf
		for c in node.child
			make_mutdict!(mutdict, mutloc, c, gaps)
		end
	end
end

"""
	find_mut_root(n::TreeTools.TreeNode, i::Int64, val)

Find most recent ancestor of `n` that bears the mutation `i`-->`val`. Return its label.
"""
function find_mut_root(n::TreeTools.TreeNode, i::Int64, val)
	if n.isroot
		return n.label
	end
	for m in  n.data.mutations
		if m.i == i && m.new == val
			return n.label
		end
	end
	#
	return find_mut_root(n.anc, i, val)
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
#
# Introspecting the database.
#

"""
    introspect(conn::LibPQ.Connection) :: PGCatalog

Introspects a Postgres database and generates a model of the database
structure.
"""
function introspect(conn)
    # Create the catalog object.
    res = execute(conn, """
        SELECT current_database()
    """)
    # current_database = res[1,1]
    current_database = first(res).current_database
    cat = PGCatalog(current_database)

    # Extract schemas.
    oid2schema = Dict{UInt32,PGSchema}()
    res = execute(conn, """
        SELECT n.oid, n.nspname
        FROM pg_catalog.pg_namespace n
        ORDER BY n.nspname
    """)
    foreach(columntable(res)...) do oid, nspname
        scm = add_schema!(cat, nspname)
        oid2schema[oid] = scm
    end

    # Extract ENUM labels.
    oid2labels = Dict{UInt32,Vector{String}}()
    res = execute(conn, """
        SELECT e.enumtypid, e.enumlabel
        FROM pg_catalog.pg_enum e
        ORDER BY e.enumtypid, e.enumsortorder, e.oid
    """)
    foreach(columntable(res)...) do enumtypid, enumlabel
        lbls = get!(() -> String[], oid2labels, enumtypid)
        push!(lbls, enumlabel)
    end

    # Extract data types.
    oid2type = Dict{UInt32,PGType}()
    res = execute(conn, """
        SELECT t.oid, t.typnamespace, t.typname, t.typtype
        FROM pg_catalog.pg_type t
        ORDER BY t.typnamespace, t.typname
    """)
    foreach(columntable(res)...) do oid, typnamespace, typname, typtype
        scm = oid2schema[typnamespace]
        typ =
            if typtype == 'e'
                add_type!(scm, typname, oid2labels[oid])
            else
                add_type!(scm, typname)
            end
        oid2type[oid] = typ
    end

    # Extract tables.
    oid2table = Dict{UInt32,PGTable}()
    res = execute(conn, """
        SELECT c.oid, c.relnamespace, c.relname
        FROM pg_catalog.pg_class c
        WHERE c.relkind IN ('r', 'v') AND
              HAS_TABLE_PRIVILEGE(c.oid, 'SELECT')
        ORDER BY c.relnamespace, c.relname
    """)
    foreach(columntable(res)...) do oid, relnamespace, relname
        scm = oid2schema[relnamespace]
        tbl = add_table!(scm, relname)
        oid2table[oid] = tbl
    end

    # Extract columns.
    oidnum2column = Dict{Tuple{UInt32,Int16},PGColumn}()
    res = execute(conn, """
        SELECT a.attrelid, a.attnum, a.attname, a.atttypid, a.attnotnull
        FROM pg_catalog.pg_attribute a
        WHERE a.attnum > 0 AND
              NOT a.attisdropped
        ORDER BY a.attrelid, a.attnum
    """)
    foreach(columntable(res)...) do attrelid, attnum, attname, atttypid, attnotnull
        attrelid in keys(oid2table) || return
        tbl = oid2table[attrelid]
        typ = oid2type[atttypid]
        col = add_column!(tbl, attname, typ, attnotnull)
        oidnum2column[(attrelid, attnum)] = col
    end

    # Extract default values.
    res = execute(conn, """
        SELECT a.adrelid, a.adnum, pg_get_expr(a.adbin, a.adrelid)
        FROM pg_catalog.pg_attrdef a
        ORDER BY a.adrelid, a.adnum
    """)
    foreach(columntable(res)...) do adrelid, adnum, adsrc
        (adrelid, adnum) in keys(oidnum2column) || return
        col = oidnum2column[(adrelid, adnum)]
        set_default!(col, adsrc)
    end

    # Extract unique keys.
    oid2unique_key = Dict{UInt32,PGUniqueKey}()
    res = execute(conn, """
        SELECT c.oid, c.conname, c.contype, c.conrelid, c.conkey
        FROM pg_catalog.pg_constraint c
        WHERE c.contype IN ('p', 'u')
        ORDER BY c.conrelid, c.conname
    """)
    foreach(columntable(res)...) do oid, conname, contype, conrelid, conkey
        conrelid in keys(oid2table) || return
        tbl = oid2table[conrelid]
        cols = PGColumn[]
        for num in conkey
            num > 0 || return
            push!(cols, oidnum2column[(conrelid, num)])
        end
        primary = contype == 'p'
        uk = add_unique_key!(tbl, conname, cols, primary)
        oid2unique_key[oid] = uk
    end

    # Extract foreign keys.
    oid2foreign_key = Dict{UInt32,PGForeignKey}()
    res = execute(conn, """
        SELECT c.oid, c.conname, c.conrelid, c.conkey, c.confrelid, c.confkey,
               c.confdeltype, c.confupdtype
        FROM pg_catalog.pg_constraint c
        WHERE c.contype = 'f'
        ORDER BY c.conrelid, c.conname
    """)
    foreach(columntable(res)...) do oid, conname, conrelid, conkey, confrelid, confkey, confdeltype, confupdtype
        conrelid in keys(oid2table) || return
        tbl = oid2table[conrelid]
        cols = PGColumn[]
        for num in conkey
            num > 0 || return
            push!(cols, oidnum2column[(conrelid, num)])
        end
        confrelid in keys(oid2table) || return
        ttbl = oid2table[confrelid]
        tcols = PGColumn[]
        for num in confkey
            num > 0 || return
            push!(tcols, oidnum2column[(confrelid, num)])
        end
        on_delete = confdeltype == 'a' ? "NO ACTION" :
                    confdeltype == 'r' ? "RESTRICT" :
                    confdeltype == 'c' ? "CASCADE" :
                    confdeltype == 'n' ? "SET NULL" :
                    confdeltype == 'd' ? "SET DEFAULT" : ""
        on_update = confupdtype == 'a' ? "NO ACTION" :
                    confupdtype == 'r' ? "RESTRICT" :
                    confupdtype == 'c' ? "CASCADE" :
                    confupdtype == 'n' ? "SET NULL" :
                    confupdtype == 'd' ? "SET DEFAULT" : ""
        fk = add_foreign_key!(tbl, conname, cols, ttbl, tcols, on_delete, on_update)
        oid2foreign_key[oid] = fk
    end

    # Extract comments
    res = execute(conn, """
        SELECT c.relname, d.objoid, d.objsubid, d.description
        FROM pg_catalog.pg_description d JOIN
             pg_catalog.pg_class c ON (d.classoid = c.oid)
        WHERE d.classoid IN ('pg_catalog.pg_namespace'::regclass,
                             'pg_catalog.pg_type'::regclass,
                             'pg_catalog.pg_class'::regclass,
                             'pg_catalog.pg_constraint'::regclass)
        ORDER BY d.objoid, d.classoid, d.objsubid
    """)
    foreach(columntable(res)...) do relname, objoid, objsubid, description
        if relname == "pg_namespace"
            scm = oid2schema[objoid]
            set_comment!(scm, description)
        elseif relname == "pg_type"
            typ = oid2type[objoid]
            set_comment!(typ, description)
        elseif relname == "pg_class" && (objoid, objsubid) in keys(oidnum2column)
            col = oidnum2column[(objoid, objsubid)]
            set_comment!(col, description)
        elseif relname == "pg_class" && objsubid == 0 && objoid in keys(oid2table)
            tbl = oid2table[objoid]
            set_comment!(tbl, description)
        elseif relname == "pg_constraint" && objoid in keys(oid2unique_key)
            uk = oid2unique_key[objoid]
            set_comment!(uk, description)
        elseif relname == "pg_constraint" && objoid in keys(oid2foreign_key)
            fk = oid2foreign_key[objoid]
            set_comment!(fk, description)
        end
    end

    cat
end


text:
####
#### Coverage summary, printed as "(percentage) covered".
####
#### Useful for CI environments that just want a summary (eg a Gitlab setup).
####

using Coverage
coverage = cd(joinpath(@__DIR__, "..", "..")) do
    coverage = process_folder()
    covered_lines, total_lines = get_summary(coverage)
    percentage = covered_lines / total_lines * 100
    println("($(percentage)%) covered")

    coverage
end
LCOV.writefile("coverage-lcov.info", coverage)

cd(joinpath(@__DIR__, "..", "..")) do
    run(`genhtml test/coverage/coverage-lcov.info --output-directory test/coverage/html`)
end