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

using Markdown
using InteractiveUtils

# ╔═╡ 77311686-7a2b-11eb-2bc7-95b2a2211969
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

# ╔═╡ 0d124834-7ef3-11eb-1d64-d7e7ff02159c
begin
	const DF = AbstractDataFrame # Union{DataFrame, SubDataFrame}

	include("./utils.jl");
	include("./dt_utils.jl");
end

# ╔═╡ 5435ff66-7a2b-11eb-1474-b96dcad21315
md"""
## C03w04: Decision Trees in Practice

In this assignment we will explore various techniques for preventing overfitting in decision trees. We will extend the implementation of the binary decision trees that we implemented in the previous assignment. <br/>
We will have use our solutions from the previous assignment and extend them.

In this assignment we will:
  -  Implement binary decision trees with different early stopping methods.
  -  Compare models with different stopping parameters.
  -  Visualize the concept of overfitting in decision trees.
"""

# ╔═╡ 771130dc-7a2b-11eb-170f-934c11fceede
md"""
#### Load LendingClub dataset
"""

# ╔═╡ 76f62e40-7a2b-11eb-095e-f36e207f06a2
begin
	loans =  CSV.File(("../../ML_UW_Spec/C03/data/lending-club-data.csv");
	header=true) |> DataFrame;
	size(loans)
end

# ╔═╡ 76de16d4-7a2b-11eb-09c6-e36b6c893c1f
md"""
Like the previous assignment, we reassign the labels to have +1 for a safe loan, and -1 for a risky (bad) loan.
"""

# ╔═╡ 76c17c4a-7a2b-11eb-2c47-97106bd58f99
insertcols!(loans, :safe_loans => ifelse.(loans.bad_loans .== 0, 1, -1),
 	makeunique=true);

# ╔═╡ 76a83684-7a2b-11eb-1961-cb133db8c164
md"""
Unlike the previous assignment where we used several features, in this assignment, we will just be using 4 categorical
features:

  1. grade of the loan
  2. the length of the loan term
  3. the home ownership status: own, mortgage, rent
  4. number of years of employment.

Since we are building a binary decision tree, we will have to convert these categorical features to a binary representation in a subsequent section using 1-hot encoding.
"""

# ╔═╡ 768d6b08-7a2b-11eb-388a-8763c8da0a51
begin
  	const Features = [
		:grade,          # grade of the loan
		:emp_length,     # number of years of employment
		:home_ownership, # home_ownership status: own, mortgage or rent
		:term,           # the term of the loan
	]

	const Target = :safe_loans # prediction target (y) (+1 means safe, -1 is risky)

	# Extract the feature columns and target column
	select!(loans, [Features..., Target]);
	length(names(loans)), names(loans)
end

# ╔═╡ 7673d8e4-7a2b-11eb-02b4-2f837181b4e9
md"""
#### Subsample dataset to make sure classes are balanced

Just as we did in the previous assignment, we will undersample the larger class (safe loans) in order to balance out our dataset. This means we are throwing away many data points.
"""

# ╔═╡ 1f459740-7ef3-11eb-34ff-9b3613ca036c
begin
	len_loans₀	 = size(loans)[1]
  	safe_loans₀, risky_loans₀ = df_partition(loans, Target)
  	p_safe_loans₀, p_risky_loans₀ = df_ratios(safe_loans₀, risky_loans₀, len_loans₀);

  	with_terminal() do
    	@printf("Number of safe loans  : %d\n", size(safe_loans₀)[1])
    	@printf("Number of risky loans : %d\n", size(risky_loans₀)[1])
    	@printf("Percentage of safe loans:  %1.2f%%\n", p_safe_loans₀ * 100.)
    	@printf("Percentage of risky loans: %3.2f%%\n", p_risky_loans₀ * 100.)
  	end
end

# ╔═╡ 32a93648-7ef3-11eb-071c-0fb7dd06e885
loans_data = resample!(safe_loans₀, risky_loans₀);

# ╔═╡ 31f2811e-7ef3-11eb-2890-7577dda5d4de
begin
	len_loans₁	 = size(loans_data)[1]
  	safe_loans₁, risky_loans₁ = df_partition(loans_data, Target)
  	p_safe_loans₁, p_risky_loans₁ = df_ratios(safe_loans₁, risky_loans₁, len_loans₁);

 	with_terminal() do
    	@printf("Number of safe loans  : %d\n", size(safe_loans₁)[1])
    	@printf("Number of risky loans : %d\n", size(risky_loans₁)[1])
    	@printf("Percentage of safe loans:  %1.2f%%\n", p_safe_loans₁ * 100.)
    	@printf("Percentage of risky loans: %3.2f%%\n", p_risky_loans₁ * 100.)
  	end
end

# ╔═╡ 760c171a-7a2b-11eb-1341-19ba5a1824a1
md"""
#### Transform categorical data into binary features

In this assignment, we will implement **binary decision trees** (decision trees for binary features, a specific case of categorical variables taking on two values, e.g., true/false). Since all of our features are currently categorical features, we want to turn them into binary features.

For instance, the `home_ownership` feature represents the home ownership status of the loanee, which is either `own`, `mortgage` or `rent`. For example, if a data point has the feature
```
   {'home_ownership': 'RENT'}
```
we want to turn this into three features:
```
 {
   'home_ownership = OWN'      : 0,
   'home_ownership = MORTGAGE' : 0,
   'home_ownership = RENT'     : 1
 }
```
"""

# ╔═╡ 75d98df2-7a2b-11eb-296e-ef4d312ffcfa
loans_data[1:5, :]

# ╔═╡ 75c1a64e-7a2b-11eb-21f6-55dd3a8ce164
begin
	hot_encode!(loans_data)
	loans_data[1:5, :]
end

# ╔═╡ 75a52c44-7a2b-11eb-2f02-2312f3a3dac9
begin
	features = names(loans_data) |> a_ -> Symbol.(a_);
	deleteat!(features,
		findall(x -> x == Target, features))
end

# ╔═╡ 75736a9c-7a2b-11eb-2278-ab37447cf22d
## Expect 25 features after hot-encoding
@test length(features) == 25

# ╔═╡ 116eb5f8-7a33-11eb-2c4c-5f937540f987
length(loans_data.grade_A)
# vs 46508 in original notebook

# ╔═╡ a5c56d68-7a32-11eb-2ef6-e58ed65dc561
md"""
**Checkpoint:** Make sure the following answers match up.
"""

# ╔═╡ a5a1cdf6-7a32-11eb-18dd-0dfee8eb885c
@test sum(loans_data.grade_A) == 6537
# vs 6422 in original notebook

# ╔═╡ a56e44e0-7a32-11eb-2587-6f5f8b8c21f3
md"""
#### Train-test split

We split the data into a train test split with 80% of the data in the training set and 20% of the data in the test set.
"""

# ╔═╡ 6fad8afe-7a33-11eb-3208-6b289c3ccd23
begin
	train_data, validation_data = train_test_split(loans_data; split=0.8, seed=42);
	size(train_data), size(validation_data)
    ## vs ((37224, 26), (9284, 26)) ...
end

# ╔═╡ 6f88e4f6-7a33-11eb-10cf-e3e8608e616b
md"""
### Early stopping methods for decision trees

In this section, we will extend the **binary tree implementation** from the previous assignment in order to handle some early stopping conditions. Recall the 3 early stopping methods that were discussed in lecture:

  1. Reached a **maximum depth**. (set by parameter `max_depth`).
  2. Reached a **minimum node size**. (set by parameter `min_node_size`).
  3. Don't split if the **gain in error reduction** is too small. (set by parameter `min_error_reduction`).

For the rest of this assignment, we will refer to these three as **early stopping conditions 1, 2, and 3**.

#### Early stopping condition 1: Maximum depth

Recall that we already implemented the maximum depth stopping condition in the previous assignment. In this assignment, we will experiment with this condition a bit more and also write code to implement the 2nd and 3rd early stopping conditions.

We will be reusing code from the previous assignment and then building upon this.
"""

# ╔═╡ 2e859492-7b2b-11eb-1d64-d7e7ff02159c
md"""
#### Early stopping condition 2: Minimum node size

The function `reached_minimum_node_size` takes 2 arguments:

1. The `data` (from a node)
2. The minimum number of data points that a node is allowed to split on, `min_node_size`.

This function simply calculates whether the number of data points at a given node is less than or equal to the specified minimum node size. This function will be used to detect this early stopping condition in the `decision_tree_create` function.
"""

# ╔═╡ 49f07e5c-7b2b-11eb-34ff-9b3613ca036c
function reached_minimum_node_size(data::DF, min_node_size)
    """true when number of data points is ≤ to the minimum node size."""
    size(data)[1] ≤ min_node_size
end

# ╔═╡ 84eab9ca-7b2b-11eb-2890-7577dda5d4de
md"""
**Quiz Question:** Given an intermediate node with 6 safe loans and 3 risky loans, if the `min_node_size` parameter is 10, what should the tree learning algorithm do next?
  - it should stop (early), creating a leaf and return it
"""

# ╔═╡ 94fccdee-7b2b-11eb-071c-0fb7dd06e885
md"""
#### Early stopping condition 3: Minimum gain in error reduction

The function `error_reduction` takes 2 arguments:

  1. The error **before** a split, `error_before_split`.
  2. The error **after** a split, `error_after_split`.

This function computes the gain in error reduction, i.e., the difference between the error before the split and that after the split. This function will be used to detect this early stopping condition in the `decision_tree_create` function.
"""

# ╔═╡ a58dfb04-7b2b-11eb-30dd-677a45cd9920
function error_reduction(error_before_split, error_after_split)
	"""Return the error before the split minus the error after the split."""
	error_before_split - error_after_split
end

# ╔═╡ c4c23a3c-7b2b-11eb-37bc-7df8411fb645
md"""
**Quiz Question:** Assume an intermediate node has 6 safe loans and 3 risky loans.  For each of 4 possible features to split on, the error reduction is 0.0, 0.05, 0.1, and 0.14, respectively. If the **minimum gain in error reduction** parameter is set to 0.2, what should the tree learning algorithm do next?
  - it should stop (early), creating a leaf and return it
"""

# ╔═╡ ce3be5b8-7b2b-11eb-1bbb-e39f6ba35768
md"""
#### Grabbing binary decision tree helper functions from past assignment

Recall from the previous assignment that we wrote a function `inter_node_num_mistakes` that calculates the number of **misclassified examples** when predicting the **majority class**. This is used to help determine which feature is best to split on at a given node of the tree.

cf. include
"""

# ╔═╡ 6f3011da-7a33-11eb-2cd0-b7bb04b7b545
## TODO move...

with_terminal() do
	@testset "begin sanity checks" begin
		## Test case 1
		ex1_labels = DataFrame(:x => [-1, -1, 1, 1, 1])
		@test inter_node_num_mistakes(ex1_labels.x) == 2

		## Test case 2
		ex2_labels = [-1, -1, 1, 1, -1, 1, 1]
		@test inter_node_num_mistakes(ex2_labels) == 3

		## Test case 3
		ex3_labels = DataFrame(:x => [-1, -1, -1, -1, 1, -1, 1])
		@test inter_node_num_mistakes(ex3_labels.x) == 2
	end
end

# ╔═╡ ce1d2376-7b2b-11eb-14c6-3dd9dd0d2303
md"""
We then wrote a function `best_splitting_feature` that finds the best feature to split on given the data and a list of features to consider.

cf. include
"""

# ╔═╡ ce025442-7b2b-11eb-31b8-69bdfe81ed69
@test best_splitting_feature(train_data, features, Target) == Symbol("term_ 60 months")

# ╔═╡ cde76c36-7b2b-11eb-2ad1-017f6ac940dd
md"""
Finally, recall the function create_leaf from the previous assignment, which creates a leaf node given a set of target values.

cf. include
"""

# ╔═╡ cdcdcc7c-7b2b-11eb-3b98-bd0e15b4de11
md"""
### Incorporating new early stopping conditions in binary decision tree implementation

Now, you will implement a function that builds a decision tree handling the three early stopping conditions described in this assignment.  In particular, you will write code to detect early stopping conditions 2 and 3.  You implemented above the functions needed to detect these conditions.  The 1st early stopping condition, **max_depth**, was implemented in the previous assigment and you will not need to reimplement this.  In addition to these early stopping conditions, the typical stopping conditions of having no mistakes or no more features to split on (which we denote by "stopping conditions" 1 and 2) are also included as in the previous assignment.

**Implementing early stopping condition 2: minimum node size:**

* **Step 1:** Use the function **reached_minimum_node_size** that you implemented earlier to write an if condition to detect whether we have hit the base case, i.e., the node does not have enough data points and should be turned into a leaf. Don't forget to use the `min_node_size` argument.
* **Step 2:** Return a leaf. This line of code should be the same as the other (pre-implemented) stopping conditions.

**Implementing early stopping condition 3: minimum error reduction:**

**Note:** This has to come after finding the best splitting feature so we can calculate the error after splitting in order to calculate the error reduction.

* **Step 1:** Calculate the **classification error before splitting**.  Recall that classification error is defined as:

$$\text{classification error} = \frac{\text{\# mistakes}}{\text{\# total examples}}$$

* **Step 2:** Calculate the **classification error after splitting**. This requires calculating the number of mistakes in the left and right splits, and then dividing by the total number of examples.
* **Step 3:** Use the function **error_reduction** to that you implemented earlier to write an if condition to detect whether  the reduction in error is less than the constant provided (`min_error_reduction`). Don't forget to use that argument.
* **Step 4:** Return a leaf. This line of code should be the same as the other (pre-implemented) stopping conditions.

"""

# ╔═╡ b6fc15ae-7a37-11eb-071c-0fb7dd06e885
function decision_tree_create(df::DF, features, target;
		curr_depth=0, max_depth=10, min_node_size=1, min_error_reduction=0.0,
		verbose=false) where {DF <: Union{DataFrame, SubDataFrame}}

    rem_features = copy(features) # Make a copy of the features.

    target_vals = df[!, target]
    if verbose
      println("-------------------------------------------------------------")
      @printf("Subtree, depth = %s (%s data points).\n", curr_depth,
			size(target_vals)[1])
	end

    if inter_node_num_mistakes(target_vals) == 0
		## Stopping cond1: no mistakes at current node?
        verbose && @printf("Stopping condition 1 reached.")
        return create_leaf(target_vals)

	elseif length(rem_features) == 0
    	## Stopping cond2: no remaining features to consider splitting on?
        verbose && @printf("Stopping condition 2 reached.")
        return create_leaf(target_vals)

	elseif curr_depth ≥ max_depth
    	## Additional stopping condition (limit tree depth)
        verbose && @printf("Early stopping condition 1 reached. Maximum depth.")
        return create_leaf(target_vals)

	elseif reached_minimum_node_size(df, min_node_size)
		## Early stopping condition 2:
        verbose && @printf("Early stopping condition 2 reached. Minimum node size.")
        return create_leaf(target_vals)
	end

    ## Now Find the best splitting feature
    splitting_feature = best_splitting_feature(df, features, target)

    # Split on the best feature that we found.
    l_split = df[df[!, splitting_feature] .== 0, :]
    r_split = df[df[!, splitting_feature] .== 1, :]

	## Early stopping condition 3: Minimum error reduction (prep)
	err_bef_split = inter_node_num_mistakes(target_vals) / size(df)[1]

	## Calculate the error after splitting
    l_mistakes = inter_node_num_mistakes(l_split[!, target])
    r_mistakes = inter_node_num_mistakes(r_split[!, target])
    err_aft_split = (l_mistakes + r_mistakes) / size(df)[1]

	if error_reduction(err_bef_split, err_aft_split) < min_error_reduction
        verbose && @printf("Early stopping condition 3. Minimum error reduction.")
        return create_leaf(target_vals)
	end

	deleteat!(rem_features,
		findall(x -> x == splitting_feature, rem_features))

    verbose && @printf("Split on feature %s. (%s, %s)\n", splitting_feature,
		size(l_split)[1], size(r_split)[1])

    # Create a leaf node if the split is "perfect"
    if size(l_split)[1] == size(df)[1]
        verbose && println("Creating leaf node.")
        return create_leaf(l_split[!, target])
	end

    if size(r_split)[1] == size(df)[1]
        verbose && println("Creating leaf node.")
        return create_leaf(r_split[!, target])
	end

    # Recurse on left and right subtrees
    ltree = decision_tree_create(l_split, rem_features, target;
						curr_depth=curr_depth + 1, max_depth,
						min_node_size, min_error_reduction, verbose)
    rtree = decision_tree_create(r_split, rem_features, target;
						curr_depth=curr_depth + 1, max_depth,
						min_node_size, min_error_reduction, verbose)

    return create_leaf(target_vals;
						splitting_feature=string(splitting_feature),
						left=ltree,
						right=rtree,
						is_leaf=false)
end

# ╔═╡ cdb69d9c-7b2b-11eb-12fa-4d42c67876b2
md"""

Run the following test code to check the implementation.
"""

# ╔═╡ cd9e2b4a-7b2b-11eb-20d9-87257af20d51
begin
	small_data_decision_tree = decision_tree_create(train_data, features, :safe_loans;
		max_depth=2, min_node_size=10, min_error_reduction=0.0, verbose=true)

	@test count_nodes(small_data_decision_tree) == 7
end

# ╔═╡ 753b3df2-7a2b-11eb-28c0-b5bdffef414b
md"""
### Build the tree!

Let's now train a tree model **ignoring early stopping conditions 2 and 3** so that we get the same tree as in the previous assignment.

To ignore these conditions, we set `min_node_size=0` and `min_error_reduction=-1` (a negative value).

Call this tree `dt6_model`.

"""

# ╔═╡ be0730f6-7a3d-11eb-30dd-677a45cd9920
## first no early stopping for later comparison
dt6_model_nes = decision_tree_create(train_data, features, :safe_loans;
	max_depth=6, min_node_size=0, min_error_reduction=-1.)

# ╔═╡ a92a0af0-7b34-11eb-2890-7577dda5d4de
md"""
Now that your code is working, we will train a tree model on the **train_data** with
  - `max_depth = 6`
  - `min_node_size = 100`,
  - `min_error_reduction = 0.0`
"""

# ╔═╡ c88e3ac2-7b34-11eb-0a69-8fbbd0859786
dt6_model = decision_tree_create(train_data, features, :safe_loans;
	max_depth=6, min_node_size=100, min_error_reduction=0.0)

# ╔═╡ c872733c-7b34-11eb-05df-8b16a1439d68
md"""
#### Making Predictions
"""

# ╔═╡ c85ab936-7b34-11eb-37bc-7df8411fb645
md"""
Recall that in the previous assignment you implemented a function `classify` to classify a new point `x` using a given `tree`.

cf. include

Now, let's consider the first example of the validation set and see what the dt6_model model predicts for this data point.
"""

# ╔═╡ c8405834-7b34-11eb-30dd-677a45cd9920
first_dp = validation_data[1, :]

# ╔═╡ 594b159e-7b35-11eb-071c-0fb7dd06e885
with_terminal() do
	println("Predicted class: $(classify(dt6_model, first_dp))")
end

# ╔═╡ 59304232-7b35-11eb-2890-7577dda5d4de
with_terminal() do
	classify(dt6_model, first_dp, annotate=true)
end

# ╔═╡ 59141350-7b35-11eb-34ff-9b3613ca036c
md"""
Let's now recall the prediction path for the decision tree learned in the previous assignment, which we recreated here as dt6_model_nes
"""

# ╔═╡ aca4b40c-7b35-11eb-20d9-87257af20d51
with_terminal() do
	classify(dt6_model_nes, first_dp, annotate=true)
end

# ╔═╡ ac83767a-7b35-11eb-0a51-851d7325ad47
md"""
**Quiz Question:** For `dt6_model` trained with `max_depth=6`, `min_node_size=100`, `min_error_reduction=0.0`, is the prediction path for `validation_set[0]` shorter, longer, or the same as for `dt6_model_nes` that ignored the early stopping conditions 2 and 3?
  - it is the same

**Quiz Question:** For `dt6_model` trained with `max_depth=6`, `min_node_size=100`, `min_error_reduction=0.0`, is the prediction path for **any point** always shorter, always longer, always the same, shorter or the same, or longer or the same as for `dt6_model_nes` that ignored the early stopping conditions 2 and 3?
  - It should be shorter or the same

**Quiz Question:** For a tree trained on **any** dataset using `max_depth=6`, `min_node_size=100`, `min_error_reduction=0.0`, what is the maximum number of splits encountered while making a single prediction?
  - max-depth, hence 6
"""

# ╔═╡ ac686542-7b35-11eb-08e3-635d3737c018
md"""
#### Evaluating the model

Now let us evaluate the model that we have trained. You implemented this evaluation in the function `eval_classification_error` from the previous assignment.

Now, let's use this function to evaluate the classification error of `dt6_model` on the **validation_data**.
"""

# ╔═╡ ac4e3b42-7b35-11eb-0a69-8fbbd0859786
with_terminal() do
	@printf("%1.4f\n",
		round(eval_classification_error(dt6_model, validation_data, Target),
			digits=4))
end

# ╔═╡ c825a994-7b34-11eb-071c-0fb7dd06e885
md"""
Now, evaluate the validation error using `dt6_model_nes`.
"""

# ╔═╡ 3bf28f8e-7b37-11eb-2890-7577dda5d4de
with_terminal() do
	@printf("%1.4f\n",
		round(eval_classification_error(dt6_model_nes, validation_data, Target),
			digits=4))
end

# ╔═╡ 58baafe8-7b37-11eb-071c-0fb7dd06e885
md"""
**Quiz Question:** Is the validation error of the new decision tree (using early stopping conditions 2 and 3) lower than, higher than, or the same as that of the old decision tree from the previous assignment?
  - should be lower! (pretty close)...
"""

# ╔═╡ 7bb705f0-7b37-11eb-0a51-851d7325ad47
md"""
### Exploring the effect of max_depth

We will compare three models trained with different values of the stopping criterion. We intentionally picked models at the extreme ends (**too small**, **just right**, and **too large**).

Train three models with these parameters:

1. **model_1**: max_depth = 2 (too small)
2. **model_2**: max_depth = 6 (just right)
3. **model_3**: max_depth = 14 (may be too large)

For each of these three, we set `min_node_size = 0` and `min_error_reduction = -1`.

**Note:** Each tree can take up to a few minutes to train. In particular, `model_3` will probably take the longest to train.

"""

# ╔═╡ 7ba0b250-7b37-11eb-08e3-635d3737c018
begin
	model_1 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=2, min_node_size=0, min_error_reduction=-1)
	model_2 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=6, min_node_size=0, min_error_reduction=-1)
	model_3 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=14, min_node_size=0, min_error_reduction=-1)
end

# ╔═╡ 7b815040-7b37-11eb-0a69-8fbbd0859786
md"""
#### Evaluating the models

Let us evaluate the models on the train and validation data. Let us start by evaluating the classification error on the training data:
"""

# ╔═╡ 7b68e74e-7b37-11eb-05df-8b16a1439d68
function print_eval(models, df, target;
		df_label=:training, offset=0, digits=4)
	for (ix, m) ∈ enumerate(models)
		cr = eval_classification_error(m, df, target)
		@printf("%10s data, classification error (model %2d): %1.4f\n",
			df_label, ix + offset, round(cr; digits=digits))
	end
end

# ╔═╡ 7b54fec0-7b37-11eb-37bc-7df8411fb645
with_terminal() do
	print_eval([model_1, model_2, model_3], train_data, Target)
end

# ╔═╡ 7b25f5d8-7b37-11eb-30dd-677a45cd9920
md"""
Now evaluate the classification error on the validation data.
"""

# ╔═╡ 69b44c90-7b38-11eb-14c6-3dd9dd0d2303
with_terminal() do
	print_eval([model_1, model_2, model_3], validation_data, Target)
end

# ╔═╡ 6993ba5c-7b38-11eb-31b8-69bdfe81ed69
md"""
**Quiz Question:** Which tree has the smallest error on the validation data?
  - model2 (in this case, with this sample)

**Quiz Question:** Does the tree with the smallest error in the training data also have the smallest error in the validation data?
  - no (in this case, with this sample)

**Quiz Question:** Is it always true that the tree with the lowest classification error on the **training** set will result in the lowest classification error in the **validation** set?
  - no

### Measuring the complexity of the tree

Recall in the lecture that we talked about deeper trees being more complex. We will measure the complexity of the tree as

```
  complexity(T) = number of leaves in the tree T
```
"""

# ╔═╡ 6974e686-7b38-11eb-2ad1-017f6ac940dd
md"""
Compute the number of nodes in model\_1, model\_2, and model\_3.
"""

# ╔═╡ cccf9118-7b38-11eb-2f79-f9575218f0a5
with_terminal() do
	for (ix, m) ∈ enumerate([model_1, model_2, model_3])
		cn = count_leaves(m)
  		@printf("for model %2d - num of node is: %4d\n", ix, cn)
	end
end

# ╔═╡ ccb33ee6-7b38-11eb-072e-f31cc84a982a
md"""
**Quiz Question:** Which tree has the largest complexity?
  - model 3

**Quiz Question:** Is it always true that the most complex tree will result in the lowest classification error in the **validation_set**?
  - no (most likely it will overfit the training set, and wil generalize poorly).
"""

# ╔═╡ cc976a84-7b38-11eb-12ec-a5cf5ea3c4c1
md"""
### Exploring the effect of min_error

We will compare three models trained with different values of the stopping criterion. We intentionally picked models at the extreme ends (**negative**, **just right**, and **too positive**).

Train three models with these parameters:
  1. **model_4**: `min_error_reduction = -1` (ignoring this early stopping condition)
  1. **model_5**: `min_error_reduction = 0` (just right)
  1. **model_6**: `min_error_reduction = 5` (too positive)

For each of these three, we set `max_depth = 6`, and `min_node_size = 0`.
"""

# ╔═╡ cc7c625e-7b38-11eb-1bbb-e39f6ba35768
begin
	model_4 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=6, min_node_size=0, min_error_reduction=-1)

	model_5 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=6,  min_node_size=0, min_error_reduction=0)

	model_6 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=6, min_node_size=0, min_error_reduction=5)
end

# ╔═╡ 695bb3dc-7b38-11eb-3b98-bd0e15b4de11
md"""
Calculate the accuracy of each model (model\_4, model\_5, or model\_6) on the validation set.
"""

# ╔═╡ fca4be14-7b39-11eb-2890-7577dda5d4de
with_terminal() do
	print_eval([model_4, model_5, model_6], validation_data, Target;
	df_label=:validation, offset=4)
end

# ╔═╡ fca323ea-7b39-11eb-34ff-9b3613ca036c
md"""
Using the count_leaves function, compute the number of leaves in each of each models in (model\_4, model\_5, and model\_6).
"""

# ╔═╡ fca2a480-7b39-11eb-1d64-d7e7ff02159c
with_terminal() do
	for (ix, m) in enumerate([model_4, model_5, model_6])
  		cn = count_leaves(m)
  		@printf("for model %2d - num of node is: %4d\n", ix, cn)
	end
end

# ╔═╡ 41f440dc-7b3a-11eb-08e3-635d3737c018
md"""
**Quiz Question:** Using the complexity definition above, which model (**model_4**, **model_5**, or **model_6**) has the largest complexity?
  - model_4 and model_5

  Did this match your expectation?
  - pretty much

**Quiz Question:** **model_4** and **model_5** have similar classification error on the validation set but **model_5** has lower complexity. Should you pick **model_5** over **model_4**?
  - yes
"""

# ╔═╡ 41d5b2ca-7b3a-11eb-0a69-8fbbd0859786
md"""
### Exploring the effect of min_node_size

We will compare three models trained with different values of the stopping criterion. Again, intentionally picked models at the extreme ends (**too small**, **just right**, and **just right**).

Train three models with these parameters:
  1. **model_7**: min_node_size = 0 (too small)
  2. **model_8**: min_node_size = 2000 (just right)
  3. **model_9**: min_node_size = 50000 (too large)

For each of these three, we set `max_depth = 6`, and `min_error_reduction = -1`.
"""

# ╔═╡ 41bd2950-7b3a-11eb-05df-8b16a1439d68
begin
	model_7 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=6, min_node_size=0, min_error_reduction=-1);

	model_8 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=6, min_node_size=2000, min_error_reduction=-1);

	model_9 = decision_tree_create(train_data, features, :safe_loans;
		max_depth=6, min_node_size=50000, min_error_reduction=-1)
end

# ╔═╡ 419fec1c-7b3a-11eb-37bc-7df8411fb645
md"""
Now, let us evaluate the models (model\_7, model\_8, or model\_9) on the validation_set.
"""

# ╔═╡ b571a75c-7b3a-11eb-3b98-bd0e15b4de11
with_terminal() do
	print_eval([model_7, model_8, model_9], validation_data, Target; 	df_label=:validation, offset=6)
end

# ╔═╡ b5553504-7b3a-11eb-12fa-4d42c67876b2
md"""
Using the count_leaves function, compute the number of leaves in each of each models (model\_7, model\_8, and model\_9).
"""

# ╔═╡ b539906a-7b3a-11eb-20d9-87257af20d51
with_terminal() do
	for (ix, m) ∈ enumerate([model_7, model_8, model_9])
  		cn = count_leaves(m)
  		@printf("for model %2d - num of node is: %4d\n", ix, cn)
	end
end

# ╔═╡ b51dd136-7b3a-11eb-0a51-851d7325ad47
md"""
**Quiz Question:** Using the results obtained in this section, which model (**model_7**, **model_8**, or **model_9**) would you choose to use?
  - model_8 (although not clear cut!)
"""

# ╔═╡ Cell order:
# ╟─5435ff66-7a2b-11eb-1474-b96dcad21315
# ╠═77311686-7a2b-11eb-2bc7-95b2a2211969
# ╠═0d124834-7ef3-11eb-1d64-d7e7ff02159c
# ╟─771130dc-7a2b-11eb-170f-934c11fceede
# ╠═76f62e40-7a2b-11eb-095e-f36e207f06a2
# ╟─76de16d4-7a2b-11eb-09c6-e36b6c893c1f
# ╠═76c17c4a-7a2b-11eb-2c47-97106bd58f99
# ╟─76a83684-7a2b-11eb-1961-cb133db8c164
# ╠═768d6b08-7a2b-11eb-388a-8763c8da0a51
# ╟─7673d8e4-7a2b-11eb-02b4-2f837181b4e9
# ╠═1f459740-7ef3-11eb-34ff-9b3613ca036c
# ╠═32a93648-7ef3-11eb-071c-0fb7dd06e885
# ╠═31f2811e-7ef3-11eb-2890-7577dda5d4de
# ╟─760c171a-7a2b-11eb-1341-19ba5a1824a1
# ╠═75d98df2-7a2b-11eb-296e-ef4d312ffcfa
# ╠═75c1a64e-7a2b-11eb-21f6-55dd3a8ce164
# ╠═75a52c44-7a2b-11eb-2f02-2312f3a3dac9
# ╠═75736a9c-7a2b-11eb-2278-ab37447cf22d
# ╠═116eb5f8-7a33-11eb-2c4c-5f937540f987
# ╟─a5c56d68-7a32-11eb-2ef6-e58ed65dc561
# ╠═a5a1cdf6-7a32-11eb-18dd-0dfee8eb885c
# ╟─a56e44e0-7a32-11eb-2587-6f5f8b8c21f3
# ╠═6fad8afe-7a33-11eb-3208-6b289c3ccd23
# ╟─6f88e4f6-7a33-11eb-10cf-e3e8608e616b
# ╟─2e859492-7b2b-11eb-1d64-d7e7ff02159c
# ╠═49f07e5c-7b2b-11eb-34ff-9b3613ca036c
# ╟─84eab9ca-7b2b-11eb-2890-7577dda5d4de
# ╟─94fccdee-7b2b-11eb-071c-0fb7dd06e885
# ╠═a58dfb04-7b2b-11eb-30dd-677a45cd9920
# ╟─c4c23a3c-7b2b-11eb-37bc-7df8411fb645
# ╟─ce3be5b8-7b2b-11eb-1bbb-e39f6ba35768
# ╠═6f3011da-7a33-11eb-2cd0-b7bb04b7b545
# ╟─ce1d2376-7b2b-11eb-14c6-3dd9dd0d2303
# ╠═ce025442-7b2b-11eb-31b8-69bdfe81ed69
# ╟─cde76c36-7b2b-11eb-2ad1-017f6ac940dd
# ╟─cdcdcc7c-7b2b-11eb-3b98-bd0e15b4de11
# ╠═b6fc15ae-7a37-11eb-071c-0fb7dd06e885
# ╟─cdb69d9c-7b2b-11eb-12fa-4d42c67876b2
# ╠═cd9e2b4a-7b2b-11eb-20d9-87257af20d51
# ╟─753b3df2-7a2b-11eb-28c0-b5bdffef414b
# ╠═be0730f6-7a3d-11eb-30dd-677a45cd9920
# ╟─a92a0af0-7b34-11eb-2890-7577dda5d4de
# ╠═c88e3ac2-7b34-11eb-0a69-8fbbd0859786
# ╟─c872733c-7b34-11eb-05df-8b16a1439d68
# ╟─c85ab936-7b34-11eb-37bc-7df8411fb645
# ╠═c8405834-7b34-11eb-30dd-677a45cd9920
# ╠═594b159e-7b35-11eb-071c-0fb7dd06e885
# ╠═59304232-7b35-11eb-2890-7577dda5d4de
# ╟─59141350-7b35-11eb-34ff-9b3613ca036c
# ╠═aca4b40c-7b35-11eb-20d9-87257af20d51
# ╟─ac83767a-7b35-11eb-0a51-851d7325ad47
# ╟─ac686542-7b35-11eb-08e3-635d3737c018
# ╠═ac4e3b42-7b35-11eb-0a69-8fbbd0859786
# ╟─c825a994-7b34-11eb-071c-0fb7dd06e885
# ╠═3bf28f8e-7b37-11eb-2890-7577dda5d4de
# ╟─58baafe8-7b37-11eb-071c-0fb7dd06e885
# ╟─7bb705f0-7b37-11eb-0a51-851d7325ad47
# ╠═7ba0b250-7b37-11eb-08e3-635d3737c018
# ╟─7b815040-7b37-11eb-0a69-8fbbd0859786
# ╠═7b68e74e-7b37-11eb-05df-8b16a1439d68
# ╠═7b54fec0-7b37-11eb-37bc-7df8411fb645
# ╟─7b25f5d8-7b37-11eb-30dd-677a45cd9920
# ╠═69b44c90-7b38-11eb-14c6-3dd9dd0d2303
# ╟─6993ba5c-7b38-11eb-31b8-69bdfe81ed69
# ╟─6974e686-7b38-11eb-2ad1-017f6ac940dd
# ╠═cccf9118-7b38-11eb-2f79-f9575218f0a5
# ╟─ccb33ee6-7b38-11eb-072e-f31cc84a982a
# ╟─cc976a84-7b38-11eb-12ec-a5cf5ea3c4c1
# ╠═cc7c625e-7b38-11eb-1bbb-e39f6ba35768
# ╟─695bb3dc-7b38-11eb-3b98-bd0e15b4de11
# ╠═fca4be14-7b39-11eb-2890-7577dda5d4de
# ╟─fca323ea-7b39-11eb-34ff-9b3613ca036c
# ╠═fca2a480-7b39-11eb-1d64-d7e7ff02159c
# ╟─41f440dc-7b3a-11eb-08e3-635d3737c018
# ╟─41d5b2ca-7b3a-11eb-0a69-8fbbd0859786
# ╠═41bd2950-7b3a-11eb-05df-8b16a1439d68
# ╟─419fec1c-7b3a-11eb-37bc-7df8411fb645
# ╠═b571a75c-7b3a-11eb-3b98-bd0e15b4de11
# ╟─b5553504-7b3a-11eb-12fa-4d42c67876b2
# ╠═b539906a-7b3a-11eb-20d9-87257af20d51
# ╟─b51dd136-7b3a-11eb-0a51-851d7325ad47


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct OutboundNatRule <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function OutboundNatRule(;id=nothing)
        o = new()
        validate_property(OutboundNatRule, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type OutboundNatRule

const _property_map_OutboundNatRule = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_OutboundNatRule = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ OutboundNatRule }) = collect(keys(_property_map_OutboundNatRule))
Swagger.property_type(::Type{ OutboundNatRule }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_OutboundNatRule[name]))}
Swagger.field_name(::Type{ OutboundNatRule }, property_name::Symbol) =  _property_map_OutboundNatRule[property_name]

function check_required(o::OutboundNatRule)
    true
end

function validate_property(::Type{ OutboundNatRule }, name::Symbol, val)
end

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FGlT_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FGlT")
JLLWrappers.@generate_main_file("FGlT", UUID("34b0a953-27c5-5560-82f5-8ee8a9e462bb"))
end  # module FGlT_jll


text:
## Solve the classical GP Regression ##
mutable struct Analytic{T<:Real,Tx<:AbstractVector,Ty<:AbstractVector} <:
               AbstractInference{T}
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    nSamples::Int
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    xview::Tx
    yview::Ty
    function Analytic{T}(ϵ::T) where {T}
        return new{T,Vector{T},Vector{T}}(ϵ)
    end
    function Analytic{T}(
        ϵ::T, nSamples::Integer, xview::Tx, yview::Ty
    ) where {T,Tx<:AbstractVector,Ty<:AbstractVector}
        return new{T,Tx,Ty}(ϵ, 0, nSamples, true, xview, yview)
    end
end

"""
    Analytic(;ϵ::T=1e-5)

Analytic inference structure for solving the classical GP regression with Gaussian noise

## Keyword arguments
- `ϵ::Real` : convergence criteria (not used at the moment)
"""
Analytic

function Analytic(; ϵ::T=1e-5) where {T<:Real}
    return Analytic{T}(ϵ)
end

function Base.show(io::IO, ::Analytic)
    return print(io, "Analytic Inference")
end

function init_inference(
    i::Analytic{T}, nSamples::Integer, xview::TX, yview::TY
) where {T<:Real,TX,TY}
    return Analytic{T}(conv_crit(i), nSamples, xview, yview)
end

function analytic_updates!(m::GP{T}) where {T}
    f = getf(m)
    l = likelihood(m)
    f.post.Σ = pr_cov(f) + first(l.σ²) * I
    f.post.α .= cov(f) \ (yview(m) - pr_mean(f, xview(m)))
    if !isnothing(l.opt_noise)
        g = 0.5 * (norm(mean(f), 2) - tr(inv(cov(f))))
        Δlogσ² = Optimise.apply!(l.opt_noise, l.σ², g .* l.σ²)
        l.σ² .= exp.(log.(l.σ²) .+ Δlogσ²)
    end
end

xview(i::Analytic) = i.xview
yview(i::Analytic) = i.yview

nMinibatch(i::Analytic) = i.nSamples

getρ(::Analytic{T}) where {T} = one(T)

MBIndices(i::Analytic) = 1:nSamples(i)

isStochastic(::Analytic) = false


text:
using Compat.SparseArrays

const ZeroCones = Union{MOI.EqualTo, MOI.Zeros}
const LPCones = Union{MOI.GreaterThan, MOI.LessThan, MOI.Nonnegatives, MOI.Nonpositives}

_dim(s::MOI.AbstractScalarSet) = 1
_dim(s::MOI.AbstractVectorSet) = MOI.dimension(s)

# Computes cone dimensions
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:ZeroCones}) = ci.value
function _allocateconstraint!(cone::Cone, f, s::ZeroCones)
    ci = cone.f
    cone.f += _dim(s)
    ci
end
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:LPCones}) = ci.value
function _allocateconstraint!(cone::Cone, f, s::LPCones)
    ci = cone.l
    cone.l += _dim(s)
    ci
end
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:MOI.SecondOrderCone}) = cone.l + ci.value
function _allocateconstraint!(cone::Cone, f, s::MOI.SecondOrderCone)
    push!(cone.qa, s.dimension)
    ci = cone.q
    cone.q += _dim(s)
    ci
end
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:MOI.ExponentialCone}) = cone.l + cone.q + ci.value
function _allocateconstraint!(cone::Cone, f, s::MOI.ExponentialCone)
    ci = 3cone.ep
    cone.ep += 1
    ci
end
constroffset(instance::ECOSOptimizer, ci::CI) = constroffset(instance.cone, ci::CI)
MOIU.canallocateconstraint(::ECOSOptimizer, ::Type{<:SF}, ::Type{<:SS}) = true
function MOIU.allocateconstraint!(instance::ECOSOptimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    CI{F, S}(_allocateconstraint!(instance.cone, f, s))
end

# Build constraint matrix
scalecoef(rows, coef, minus, s) = minus ? -coef : coef
scalecoef(rows, coef, minus, s::Union{MOI.LessThan, Type{<:MOI.LessThan}, MOI.Nonpositives, Type{MOI.Nonpositives}}) = minus ? coef : -coef
_varmap(f) = map(vi -> vi.value, f.variables)
_constant(s::MOI.EqualTo) = s.value
_constant(s::MOI.GreaterThan) = s.lower
_constant(s::MOI.LessThan) = s.upper
constrrows(::MOI.AbstractScalarSet) = 1
constrrows(s::MOI.AbstractVectorSet) = 1:MOI.dimension(s)
constrrows(instance::ECOSOptimizer, ci::CI{<:MOI.AbstractScalarFunction, <:MOI.AbstractScalarSet}) = 1
constrrows(instance::ECOSOptimizer, ci::CI{<:MOI.AbstractVectorFunction, MOI.Zeros}) = 1:instance.cone.eqnrows[constroffset(instance, ci)]
constrrows(instance::ECOSOptimizer, ci::CI{<:MOI.AbstractVectorFunction, <:MOI.AbstractVectorSet}) = 1:instance.cone.ineqnrows[constroffset(instance, ci)]
matrix(data::Data, s::ZeroCones) = data.b, data.IA, data.JA, data.VA
matrix(data::Data, s::Union{LPCones, MOI.SecondOrderCone, MOI.ExponentialCone}) = data.h, data.IG, data.JG, data.VG
matrix(instance::ECOSOptimizer, s) = matrix(instance.data, s)
MOIU.canloadconstraint(::ECOSOptimizer, ::Type{<:SF}, ::Type{<:SS}) = true
MOIU.loadconstraint!(instance::ECOSOptimizer, ci, f::MOI.SingleVariable, s) = MOIU.loadconstraint!(instance, ci, MOI.ScalarAffineFunction{Float64}(f), s)
function MOIU.loadconstraint!(instance::ECOSOptimizer, ci, f::MOI.ScalarAffineFunction, s::MOI.AbstractScalarSet)
    a = sparsevec(_varmap(f), f.coefficients)
    # sparsevec combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(a)
    offset = constroffset(instance, ci)
    row = constrrows(s)
    i = offset + row
    # The ECOS format is b - Ax ∈ cone
    # so minus=false for b and minus=true for A
    setconstant = _constant(s)
    if s isa MOI.EqualTo
        instance.cone.eqsetconstant[offset] = setconstant
    else
        instance.cone.ineqsetconstant[offset] = setconstant
    end
    constant = f.constant - setconstant
    b, I, J, V = matrix(instance, s)
    b[i] = scalecoef(row, constant, false, s)
    append!(I, fill(i, length(a.nzind)))
    append!(J, a.nzind)
    append!(V, scalecoef(row, a.nzval, true, s))
end
MOIU.loadconstraint!(instance::ECOSOptimizer, ci, f::MOI.VectorOfVariables, s) = MOIU.loadconstraint!(instance, ci, MOI.VectorAffineFunction{Float64}(f), s)
# SCS orders differently than MOI the second and third dimension of the exponential cone
orderval(val, s) = val
function orderval(val, s::Union{MOI.ExponentialCone, Type{MOI.ExponentialCone}})
    val[[1, 3, 2]]
end
orderidx(idx, s) = idx
expmap(i) = (1, 3, 2)[i]
function orderidx(idx, s::MOI.ExponentialCone)
    expmap.(idx)
end
function MOIU.loadconstraint!(instance::ECOSOptimizer, ci, f::MOI.VectorAffineFunction, s::MOI.AbstractVectorSet)
    A = sparse(f.outputindex, _varmap(f), f.coefficients)
    # sparse combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(A)
    colval = zeros(Int, length(A.rowval))
    for col in 1:A.n
        colval[A.colptr[col]:(A.colptr[col+1]-1)] = col
    end
    @assert !any(iszero.(colval))
    offset = constroffset(instance, ci)
    rows = constrrows(s)
    if s isa MOI.Zeros
        instance.cone.eqnrows[offset] = length(rows)
    else
        instance.cone.ineqnrows[offset] = length(rows)
    end
    i = offset + rows
    # The ECOS format is b - Ax ∈ cone
    # so minus=false for b and minus=true for A
    b, I, J, V = matrix(instance, s)
    b[i] = scalecoef(rows, orderval(f.constant, s), false, s)
    append!(I, offset + orderidx(A.rowval, s))
    append!(J, colval)
    append!(V, scalecoef(A.rowval, A.nzval, true, s))
end

function MOIU.allocatevariables!(instance::ECOSOptimizer, nvars::Integer)
    instance.cone = Cone()
    VI.(1:nvars)
end

function MOIU.loadvariables!(instance::ECOSOptimizer, nvars::Integer)
    cone = instance.cone
    m = cone.l + cone.q + 3cone.ep
    IA = Int[]
    JA = Int[]
    VA = Float64[]
    b = zeros(cone.f)
    IG = Int[]
    JG = Int[]
    VG = Float64[]
    h = zeros(m)
    c = zeros(nvars)
    instance.data = Data(m, nvars, IA, JA, VA, b, IG, JG, VG, h, 0., c)
end

MOIU.canallocate(::ECOSOptimizer, ::MOI.ObjectiveSense) = true
function MOIU.allocate!(instance::ECOSOptimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    instance.maxsense = sense == MOI.MaxSense
end
MOIU.canallocate(::ECOSOptimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
function MOIU.allocate!(::ECOSOptimizer, ::MOI.ObjectiveFunction, ::MOI.ScalarAffineFunction) end

MOIU.canload(::ECOSOptimizer, ::MOI.ObjectiveSense) = true
function MOIU.load!(::ECOSOptimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense) end
MOIU.canload(::ECOSOptimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
function MOIU.load!(instance::ECOSOptimizer, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction)
    c0 = full(sparsevec(_varmap(f), f.coefficients, instance.data.n))
    instance.data.objconstant = f.constant
    instance.data.c = instance.maxsense ? -c0 : c0
end

function MOI.optimize!(instance::ECOSOptimizer)
    cone = instance.cone
    m = instance.data.m
    n = instance.data.n
    A = ECOS.ECOSMatrix(sparse(instance.data.IA, instance.data.JA, instance.data.VA, cone.f, n))
    b = instance.data.b
    G = ECOS.ECOSMatrix(sparse(instance.data.IG, instance.data.JG, instance.data.VG, m, n))
    h = instance.data.h
    objconstant = instance.data.objconstant
    c = instance.data.c
    instance.data = nothing # Allows GC to free instance.data before A is loaded to ECOS
    ecos_prob_ptr = ECOS.setup(n, m, cone.f, cone.l, length(cone.qa), cone.qa, cone.ep, G, A, c, h, b)
    ret_val = ECOS.solve(ecos_prob_ptr)
    ecos_prob = unsafe_wrap(Array, ecos_prob_ptr, 1)[1]
    primal    = unsafe_wrap(Array, ecos_prob.x, n)[:]
    dual_eq   = unsafe_wrap(Array, ecos_prob.y, cone.f)[:]
    dual_ineq = unsafe_wrap(Array, ecos_prob.z, m)[:]
    slack     = unsafe_wrap(Array, ecos_prob.s, m)[:]
    ECOS.cleanup(ecos_prob_ptr, 0)
    objval = (instance.maxsense ? -1 : 1) * dot(c, primal) + objconstant
    instance.sol = Solution(ret_val, primal, dual_eq, dual_ineq, slack, objval)
end


text:
__precompile__()

module LPWriter

# order, indices, weights
const SOS = Tuple{Int, Vector{Int}, Vector{Float64}}

const START_REG = r"^([\.0-9eE])"
const NAME_REG = r"([^a-zA-Z0-9\!\"\#\$\%\&\(\)\/\,\.\;\?\@\_\`\'\{\}\|\~])"

function verifyname(name::String)
    if length(name) > 16
        return false
    end
    m = match(START_REG, name)
    if !isa(m, Void)
        return false
    end
    m = match(NAME_REG, name)
    if !isa(m, Void)
        return false
    end
    return true
end

function correctname(name::String)
    if length(name) > 16
        warn("Name $(name) too long. Truncating.")
        return correctname(String(name[1:16]))
    end
    m = match(START_REG, name)
    if !isa(m, Void)
        warn("Name $(name) cannot start with a period, a number, e, or E. Removing from name.")
        return correctname(replace(name, START_REG, ""))
    end
    m = match(NAME_REG, name)
    if !isa(m, Void)
        warn("Name $(name) contains an illegal character. Removing from name.")
        return correctname(replace(name, NAME_REG, ""))
    end
    return name
end

include("writer.jl")
include("reader.jl")

end


text:
function _is_cell_done(cell)
    if cell.running_disabled
        return true
    else
        return !cell.queued && !cell.running
    end
end

function nothingstring(x::Union{Nothing,AbstractString})::Union{Nothing,String}
    return x isa Nothing ? x : string(x)::String
end

@enum OutputFormat begin
    documenter_output
    franklin_output
    html_output
end

"""
    BuildOptions(
        dir::AbstractString;
        write_files::Bool=true,
        previous_dir::Union{Nothing,AbstractString}=nothing,
        output_format::OutputFormat=html_output,
        use_distributed::Bool=true
    )

Arguments:

- `dir`:
    Directory in which the Pluto notebooks are stored.
- `write_files`:
    Write files to `joinpath(dir, "\$file.html")`.
- `previous_dir::Union{Nothing,AbstractString}=Nothing`:
    Use the output from the previous run as a cache to speed up running time.
    To use the cache, specify a directory `previous_dir::AbstractString` which contains HTML or Markdown files from a previous run.
    Specifically, files are expected to be at `joinpath(previous_dir, "\$file.html")`.
    The output from the previous run may be embedded in a larger HTML or Markdown file.
    This package will extract the original output from the full file contents.
    By default, caching is disabled.
-  `output_format`:
    What file to write the output to.
    By default this is `html_output::OutputFormat` meaning that the output of the HTML method is pure HTML.
    To generate Franklin or Documenter files, use respectively `franklin_output` or `documenter_output`.
    When `BuildOptions.write_files == true` and `output_format != html_output`, the output file has a ".md" extension instead of ".html".
- `use_distributed`:
    Whether to build the notebooks in different processes.
    By default, this is enabled just like in Pluto and the notebooks are build in parallel.
    The benefit of different processes is that things are more independent of each other.
    Unfortunately, the drawback is that compilation has to happen for each process.
    By setting this option to `false`, all notebooks are built sequentially in the same process which avoids recompilation.
    This is likely quicker in situations where there are few threads available such as GitHub Runners depending on the notebook contents.
    Beware that `use_distributed=false` will not work with Pluto's built-in package manager.
"""
struct BuildOptions
    dir::String
    write_files::Bool
    previous_dir::Union{Nothing,String}
    output_format::OutputFormat
    use_distributed::Bool

    function BuildOptions(
        dir::AbstractString;
        write_files::Bool=true,
        previous_dir::Union{Nothing,AbstractString}=nothing,
        output_format::OutputFormat=html_output,
        use_distributed::Bool=true
    )
        return new(
            string(dir)::String,
            write_files,
            nothingstring(previous_dir),
            output_format,
            use_distributed
        )
    end
end

"""
    _is_notebook_done(notebook::Notebook)

Return whether all cells in the `notebook` have executed.
This method is more reliable than using `notebook.executetoken` because Pluto.jl drops that lock also after installing packages.
"""
function _notebook_done(notebook::Notebook)
    cells = [last(elem) for elem in notebook.cells_dict]
    return all(_is_cell_done, cells)
end

function extract_previous_output(html::AbstractString)::String
    start_range = findfirst(BEGIN_IDENTIFIER, html)
    @assert !isnothing(start_range)
    start = first(start_range)

    stop_range = findfirst(END_IDENTIFIER, html)
    @assert !isnothing(stop_range)
    stop = last(stop_range)

    return html[start:stop]
end

"""
    Previous(state::State, html::String)

- `state`: Previous state.
- `text`: Either
    - HTML starting with "$BEGIN_IDENTIFIER" and ending with "$END_IDENTIFIER" or
    - Contents of a Franklin/Documenter Markdown file.
"""
struct Previous
    state::Union{State,Nothing}
    text::String
end

function Previous(text::String)
    state = contains(text, STATE_IDENTIFIER) ?
        extract_state(text) :
        nothing
    text = contains(text, BEGIN_IDENTIFIER) ?
        extract_previous_output(text) :
        ""
    return Previous(state, text)
end

function Previous(bopts::BuildOptions, in_file)
    prev_dir = bopts.previous_dir
    if isnothing(prev_dir)
        return Previous(nothing, "")
    end
    name, _ = splitext(in_file)
    ext = bopts.output_format == html_output ? ".html" : ".md"
    prev_path = joinpath(prev_dir, "$name$ext")
    if !isfile(prev_path)
        return Previous(nothing, "")
    end
    text = read(prev_path, String)
    return Previous(text)
end

function reuse_previous(previous::Previous, dir, in_file)::Bool
    in_path = joinpath(dir, in_file)
    curr = path2state(in_path)

    prev = previous.state
    isnothing(prev) && return false

    sha_match = prev.input_sha == curr.input_sha
    julia_match = prev.julia_version == curr.julia_version
    reuse = sha_match && julia_match
    return reuse
end

"""
Write to a ".html" or ".md" file depending on `HTMLOptions.output_format`.
The output file is always a sibling to the file at `in_path`.
"""
function _write_main_output(in_path, text, bopts::BuildOptions, hopts::HTMLOptions)
    ext = bopts.output_format == html_output ? ".html" : ".md"
    if bopts.write_files
        dir = dirname(in_path)
        in_file = basename(in_path)
        without_extension, _ = splitext(in_file)
        out_file = "$(without_extension)$(ext)"
        out_path = joinpath(dir, out_file)
        write(out_path, text)
    end
    return nothing
end

"Used when creating the page for the first time and to restore the cache."
function _wrap_franklin_output(html)
    return "~~~\n$(html)\n~~~"
end

"Add some style overrides to make things a bit prettier and more consistent with Pluto."
function _add_documenter_style(html)
    style = """
        <style>
            table {
                display: table !important;
                margin: 2rem auto !important;
                border-top: 2pt solid rgba(0,0,0,0.2);
                border-bottom: 2pt solid rgba(0,0,0,0.2);
            }

            pre, div {
                margin-top: 1.4rem !important;
                margin-bottom: 1.4rem !important;
            }
        </style>
        """
    return string(style, '\n', html)
end

"Used when creating the page for the first time and to restore the cache."
function _wrap_documenter_output(html)
    html = _add_documenter_style(html)
    return "```@raw html\n$(html)\n```"
end

function _outcome2text(session, prev::Previous, in_path::String, bopts, hopts)::String
    text = prev.text
    if bopts.output_format == franklin_output
        text = _wrap_franklin_output(text)
    end
    if bopts.output_format == documenter_output
        text = _wrap_documenter_output(text)
    end
    _write_main_output(in_path, text, bopts, hopts)
    return text
end

function _inject_script(html, script)
    l = length(END_IDENTIFIER)
    without_end = html[1:end-l]
    return string(without_end, '\n', script, '\n', END_IDENTIFIER)
end

function _outcome2text(session, nb::Notebook, in_path::String, bopts, hopts)::String
    while !_notebook_done(nb)
        sleep(0.1)
    end

    _throw_if_error(session, nb)

    # Grab output before changing binds via `_run_dynamic!`.
    # Otherwise, the outputs look wrong when opening a page for the first time.
    html = notebook2html(nb, in_path, hopts)

    if bopts.output_format == franklin_output
        html = _wrap_franklin_output(html)
    end
    if bopts.output_format == documenter_output
        html = _wrap_documenter_output(html)
    end

    _write_main_output(in_path, html, bopts, hopts)

    # The sleep avoids `AssertionError: will_run_code(notebook)`
    @async begin
        sleep(5)
        SessionActions.shutdown(session, nb; verbose=false)
    end

    return string(html)::String
end

"""
    build_notebooks(
        bopts::BuildOptions,
        files,
        hopts::HTMLOptions=HTMLOptions();
        session=ServerSession()
    ) -> Vector{String}

Build all `files` in `dir` in parallel.
"""
function build_notebooks(
        bopts::BuildOptions,
        files,
        hopts::HTMLOptions=HTMLOptions();
        session=ServerSession()
    )::Vector{String}

    dir = bopts.dir

    # Start all the notebooks in parallel with async enabled if `use_distributed`.
    X = map(files) do in_file
        in_path = joinpath(dir, in_file)::String
        @assert isfile(in_path) "Expected .jl file at $in_path"

        previous = Previous(bopts, in_file)
        if reuse_previous(previous, dir, in_file)
            @info "Using cache for Pluto notebook at $in_file"
            return previous
        else
            @info "Starting evaluation of Pluto notebook $in_file"
            if bopts.use_distributed
                run_async = true
            else
                run_async = false
                # `use_distributed` means mostly "whether to run in a new process".
                session.options.evaluation.workspace_use_distributed = false
            end
            nb = run_notebook!(in_path, session; hopts, run_async)
            return nb
        end
    end

    H = map(zip(files, X)) do (in_file, x)
        in_path = joinpath(dir, in_file)
        html = _outcome2text(session, x, in_path, bopts, hopts)
        return html
    end

    return H
end
precompile(build_notebooks, (BuildOptions, Vector{Any}, HTMLOptions))

function _is_pluto_file(path::AbstractString)::Bool
    first(eachline(string(path))) == "### A Pluto.jl notebook ###"
end

"""
    build_notebooks(
        bopts::BuildOptions,
        hopts::HTMLOptions=HTMLOptions()
    ) -> Vector{String}

Build all ".jl" files in `dir` in parallel.
"""
function build_notebooks(
        bopts::BuildOptions,
        hopts::HTMLOptions=HTMLOptions()
    )::Vector{String}
    dir = bopts.dir
    files = filter(readdir(dir)) do file
        path = joinpath(dir, file)
        endswith(file, ".jl") && _is_pluto_file(path) && !startswith(file, TMP_COPY_PREFIX)
    end
    return build_notebooks(bopts, files, hopts)
end
precompile(build_notebooks, (BuildOptions, HTMLOptions))


text:
struct CatView{N, T, NR} <: AbstractVector{T}
    arr::NR
    len::NTuple{N,Int}
end

## Constructors ##
@inline CatView(a::AbstractVector...) = CatView(tuple(a...))

@generated function CatView(arr::NR) where NR<:Tuple
    N = length(NR.parameters)
    array_types = eltype(NR)
    array_types <: AbstractVector || throw(DomainError(array_types, "CatView is only valid on AbstractVectors."))
    T = eltype(array_types)

    quote
    len = @ntuple $N (n)->length(arr[n])
    CatView{$N, $T, $NR}(arr, len)
    end
end

## size ##
Base.size(A::CatView) = (sum(A.len),)

## get index and set index ##

function Base.getindex(A::CatView, i::Int)
    i < 1 || i > length(A) && throw(BoundsError("index out of bounds."))

    a = 0
    b = A.len[1]
    @inbounds for j = 1:length(A.len)
        if i <= b
            return A.arr[j][i-a]
        else
            a = b
            b = b + A.len[j+1]
        end
    end
end

function Base.setindex!(A::CatView, val, i::Int)
    i < 1 || i > length(A) && throw(BoundsError("index out of bounds."))

    a = 0
    b = A.len[1]
    @inbounds for j = 1:length(A.len)
        if i <= b
            A.arr[j][i-a] = val
            return val
        else
            a = b
            b = b + A.len[j+1]
        end
    end
end

#Base.@propagate_inbounds
function Base.getindex(A::CatView, idx::Tuple{Integer,Integer})
    i,j = idx
    return A.arr[i][j]
end

#Base.@propagate_inbounds
function Base.setindex!(A::CatView, val, idx::Tuple{Integer,Integer})
    i,j = idx
    return setindex!(A.arr[i], val, j)
end

## Fast iteration ##
@generated function Base.eachindex(A::CatView{N,T}) where {N,T}
    quote
    @nexprs $N (n)->(i_n = zip(repeated(n,length(A.arr[n])),eachindex(A.arr[n])))
    flatten( (@ntuple $N (n)->i_n) )
    end
end

function Base.mapreduce(f, op, A::CatView)
    reduce(op, mapreduce(f, op, aa) for aa in A.arr)
end


text:
"""
    degrees(g[, T]; dir=:out)

Degree of each vertex. Return a vector which contains the degree of each vertex in graph `g`.

# Arguments

- `g`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from Graphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `T`: result element type of degree vector; default is the element type of `g` (optional).
- `dir`: direction of degree; should be `:in`, `:out`, or `:both` (optional).

# Examples

```jldoctest
julia> using GraphLaplacians

julia> m = [0 1 1; 1 0 0; 1 0 0];

julia> GraphLaplacians.degrees(m)
3-element Vector{Int64}:
 2
 1
 1
```
"""
function degrees(g::AbstractGraph, T::DataType=eltype(g); dir::Symbol=:out)
    adj = Graphs.adjacency_matrix(g, T; dir=dir)
    degrees(adj, T; dir=dir)
end

function degree_matrix(g::AbstractGraph, T::DataType=eltype(g); dir::Symbol=:out)
    adj = Graphs.adjacency_matrix(g, T; dir=dir)
    degree_matrix(adj, T; dir=dir)
end

# Graphs.laplacian_matrix(g::AbstractGraph, T::DataType=eltype(g); dir::Symbol=:out)

function normalized_laplacian(g::AbstractGraph, T::DataType=eltype(g);
                              dir::Symbol=:both, selfloop::Bool=false)
    adj = Graphs.adjacency_matrix(g, T)
    normalized_laplacian(adj, T, dir=dir, selfloop=selfloop)
end

function scaled_laplacian(g::AbstractGraph, T::DataType=eltype(g))
    adj = adjacency_matrix(g, T)
    scaled_laplacian(adj, T)
end


text:
# select the occupied orbitals assuming an insulator
function select_occupied_orbitals(basis::PlaneWaveBasis, ψ)
    model  = basis.model
    n_spin = model.n_spin_components
    @assert iszero(basis.model.temperature)
    n_bands = div(model.n_electrons, n_spin * filled_occupation(model), RoundUp)
    [ψk[:, 1:n_bands] for ψk in ψ]
end

# Packing routines used in direct_minimization and newton algorithms.
# They pack / unpack sets of ψ's (or compatible arrays, such as hamiltonian
# applies and gradients) to make them compatible to be used in algorithms
# from IterativeSolvers.
# Some care is needed here : some operators (for instance K in newton.jl)
# are real-linear but not complex-linear. To overcome this difficulty, instead of
# seeing them as operators from C^N to C^N, we see them as
# operators from R^2N to R^2N. In practice, this is done with the
# reinterpret function from julia.
# /!\ pack_ψ does not share memory while unpack_ψ does

reinterpret_real(x) = reinterpret(real(eltype(x)), x)
reinterpret_complex(x) = reinterpret(Complex{eltype(x)}, x)

function pack_ψ(ψ)
    # TODO as an optimization, do that lazily? See LazyArrays
    vcat([vec(ψk) for ψk in ψ]...)
end

function unpack_ψ(x, sizes_ψ)
    n_bands = sizes_ψ[1][2]
    lengths = prod.(sizes_ψ)
    ends = cumsum(lengths)
    # We unsafe_wrap the resulting array to avoid a complicated type for ψ.
    # The resulting array is valid as long as the original x is still in live memory.
    map(1:length(sizes_ψ)) do ik
        unsafe_wrap(Array{complex(eltype(x))},
                    pointer(@views x[ends[ik]-lengths[ik]+1:ends[ik]]),
                    sizes_ψ[ik])
    end
end


text:
# return the i-th bit
@inline bitat(chunk::UInt64, i::Int) = (chunk >>> (64 - i)) & 0x01 == 0x01
@inline bitat(::Type{T}, chunk::T, i::Int) where {T<:Unsigned} = (chunk >>> (sizeof(T) * 8 - i)) & 1 == 1

# assume 1 byte = 8 bits
@inline bitsof(::Type{T}) where {T<:Unsigned} = sizeof(T) * 8

@inline mod64(i) = i & 63

# make a bit mask
@inline lmask(typ::Type{T}, n::Int) where {T<:Unsigned} = typemax(typ) << (bitsof(typ) - n)
@inline rmask(typ::Type{T}, n::Int) where {T<:Unsigned} = typemax(typ) >> (bitsof(typ) - n)


text:
using Test
using MPIArrays

println(methods(CyclicRange))
cr1 = CyclicRange(start=1, stop=11, blocksize=2, n_split=3, split_index=1)
cr2 = CyclicRange(start=1, stop=11, blocksize=2, n_split=3, split_index=2)
cr3 = CyclicRange(start=1, stop=11, blocksize=2, n_split=3, split_index=3)

# a a b b c c a a b b c
@show cr1 cr2 cr3

@show length(cr1) length(cr2) length(cr3)
# nblock = 1, extra_block = 0
@test MPIArrays.number_block_local_min(cr1) == 1
@test collect(cr1) == [1, 2, 7, 8]
@test collect(cr2) == [3, 4, 9, 10]
@test collect(cr3) == [5, 6, 11]
@test length(cr1) == 4
@test length(cr2) == 4
@test length(cr3) == 3

for c in [cr1, cr2, cr3]
for index in 1:length(c)
    @show c[index]
end
end

println("local_index test")
for i in 1:11
    println("global=$i, ", MPIArrays.local_index(cr1, i))
end

# CyclicRangeがAbstractRangeを継承しているのでそのまま使える
println(cr1...)
println(vcat(cr1))
println(collect(cr1))
# show(stdout, cr)

# test 1
cyclic_partitions = MPIArrays.CyclicPartitioning(array_sizes=(9, 9), n_procs=(2, 2), blocksizes=(1, 1))
println(size(cyclic_partitions))
println(cyclic_partitions[1, 1])

println("CyclicPartitioning: local_index test")
show(stdout, "text/plain", map(x -> MPIArrays.local_index(cyclic_partitions, x), Base.Iterators.product(1:9, 1:9)))
println()

# test 2
cyclic_partitions = MPIArrays.CyclicPartitioning(array_sizes=(5, 7), n_procs=(2, 2), blocksizes=(2, 2))
println(size(cyclic_partitions))
println(cyclic_partitions[1, 1])

println("CyclicPartitioning: local_index test")
show(stdout, "text/plain", map(x -> MPIArrays.local_index(cyclic_partitions, x), Base.Iterators.product(1:5, 1:7)))
println()

show(stdout, "text/plain", collect(Base.Iterators.product(cr1, cr1)))
println("")
println("rank 0 local array size: ", length.(cyclic_partitions[1]))
println("rank 1 local array size: ", length.(cyclic_partitions[2]))
println("rank 2 local array size: ", length.(cyclic_partitions[3]))
println("rank 3 local array size: ", length.(cyclic_partitions[4]))


text:
##==============================================================================
## (==)
##==============================================================================
import Base.==
## @generated compare
# Reference https://github.com/JuliaLang/julia/issues/4648

#=
For now abstract `InferenceVariable`s are considered equal if they are the same type, dims, and manifolds (abels are deprecated)
If your implentation has aditional properties such as `DynPose2` with `ut::Int64` (microsecond time) or support different manifolds
implement compare if needed.
=#
==(a::InferenceVariable,b::InferenceVariable) = typeof(a) == typeof(b) && a.dims == b.dims && a.manifolds == b.manifolds

==(a::FactorOperationalMemory, b::FactorOperationalMemory) = typeof(a) == typeof(b)

==(a::FunctorInferenceType, b::FunctorInferenceType) = typeof(a) == typeof(b)

# Generate compares automatically for all in this union
const GeneratedCompareUnion = Union{MeanMaxPPE, VariableNodeData,
                              DFGVariable, DFGVariableSummary, SkeletonDFGVariable,
                              GenericFunctionNodeData,
                              DFGFactor, DFGFactorSummary, SkeletonDFGFactor}

@generated function ==(x::T, y::T) where T <: GeneratedCompareUnion
    ignored = []
    mapreduce(n -> :(x.$n == y.$n), (a,b)->:($a && $b), setdiff(fieldnames(x), ignored))
end

##==============================================================================
## Compare
##==============================================================================

function compareField(Allc, Bllc, syms)::Bool
    (!isdefined(Allc, syms) && !isdefined(Bllc, syms)) && return true
    !isdefined(Allc, syms) && return false
    !isdefined(Bllc, syms) && return false
    return eval(:($Allc.$syms == $Bllc.$syms))
end

"""
    $(SIGNATURES)

Compare the all fields of T that are not in `skip` for objects `Al` and `Bl`.

TODO > add to func_ref.md
"""
function compareFields(Al::T,
                       Bl::T;
                       show::Bool=true,
                       skip::Vector{Symbol}=Symbol[]  )::Bool where {T}
  for field in fieldnames(T)
    (field in skip) && continue
    tp = compareField(Al, Bl, field)
    show && @debug("  $tp : $field")==nothing
    !tp && return false
  end
  return true
end

function compareFields(Al::T,
                       Bl::T;
                       show::Bool=true,
                       skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Union{Number, AbstractString}}
  #
  return Al == Bl
end

function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Union{AbstractString,Symbol}}
  #
  return Al == Bl
end

function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Union{Array{<:Number}, Number}}
  #
  (length(Al) != length(Bl)) && return false
  return norm(Al - Bl) < 1e-6
end

function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Array}
  #
  (length(Al) != length(Bl)) && return false
  for i in 1:length(Al)
    !compareAll(Al[i],Bl[i], show=false) && return false
  end
  return true
end

"""
    $(SIGNATURES)

Recursively compare the all fields of T that are not in `skip` for objects `Al` and `Bl`.

TODO > add to func_ref.md
"""
function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Tuple}
  #
  length(Al) != length(Bl) && return false
  for i in 1:length(Al)
    !compareAll(Al[i], Bl[i], show=show, skip=skip) && return false
  end
  return true
end

function compareAll(Al::T,
                    Bl::T;
                    show::Bool=true,
                    skip::Vector{Symbol}=Symbol[]  )::Bool where {T <: Dict}
  #
  (length(Al) != length(Bl)) && return false
  for (id, val) in Al
    (Symbol(id) in skip) && continue
    !compareAll(val, Bl[id], show=show, skip=skip) && return false
  end
  return true
end

function compareAll(Al::T, Bl::T; show::Bool=true, skip::Vector{Symbol}=Symbol[])::Bool where T
  @debug "Comparing types $T:"
  # @debug "  Al = $Al"
  # @debug "  Bl = $Bl"
  !compareFields(Al, Bl, show=show, skip=skip) && return false
  for field in fieldnames(T)
    field in skip && continue
    @debug("  Checking field: $field")
    (!isdefined(Al, field) && !isdefined(Al, field)) && return true
    !isdefined(Al, field) && return false
    !isdefined(Bl, field) && return false
    Ad = eval(:($Al.$field))
    Bd = eval(:($Bl.$field))
    !compareAll(Ad, Bd, show=show, skip=skip) && return false
  end
  return true
end

#Compare VariableNodeData
function compare(a::VariableNodeData, b::VariableNodeData)
    a.val != b.val && @debug("val is not equal")==nothing && return false
    a.bw != b.bw && @debug("bw is not equal")==nothing && return false
    a.BayesNetOutVertIDs != b.BayesNetOutVertIDs && @debug("BayesNetOutVertIDs is not equal")==nothing && return false
    a.dimIDs != b.dimIDs && @debug("dimIDs is not equal")==nothing && return false
    a.dims != b.dims && @debug("dims is not equal")==nothing && return false
    a.eliminated != b.eliminated && @debug("eliminated is not equal")==nothing && return false
    a.BayesNetVertID != b.BayesNetVertID && @debug("BayesNetVertID is not equal")==nothing && return false
    a.separator != b.separator && @debug("separator is not equal")==nothing && return false
    a.initialized != b.initialized && @debug("initialized is not equal")==nothing && return false
    abs(a.inferdim - b.inferdim) > 1e-14 && @debug("inferdim is not equal")==nothing && return false
    a.ismargin != b.ismargin && @debug("ismargin is not equal")==nothing && return false
    a.dontmargin != b.dontmargin && @debug("dontmargin is not equal")==nothing && return false
    a.solveInProgress != b.solveInProgress && @debug("solveInProgress is not equal")==nothing && return false
    typeof(a.softtype) != typeof(b.softtype) && @debug("softtype is not equal")==nothing && return false
    return true
end

"""
    $SIGNATURES

Compare that all fields are the same in a `::FactorGraph` variable.
"""
function compareVariable(A::DFGVariable,
                         B::DFGVariable;
                         skip::Vector{Symbol}=Symbol[],
                         show::Bool=true,
                         skipsamples::Bool=true  )::Bool
  #
  skiplist = union([:attributes;:solverDataDict;:createdTimestamp;:lastUpdatedTimestamp],skip)
  TP = compareAll(A, B, skip=skiplist, show=show)
  varskiplist = skipsamples ? [:val; :bw] : Symbol[]
  skiplist = union([:softtype;],varskiplist)
  union!(skiplist, skip)
  TP = TP && compareAll(A.solverDataDict, B.solverDataDict, skip=skiplist, show=show)

  Ad = getSolverData(A)
  Bd = getSolverData(B)

  # TP = TP && compareAll(A.attributes, B.attributes, skip=[:softtype;], show=show)
  varskiplist = union(varskiplist, [:softtype])
  union!(varskiplist, skip)
  TP = TP && compareAll(Ad, Bd, skip=varskiplist, show=show)
  TP = TP && typeof(Ad.softtype) == typeof(Bd.softtype)
  TP = TP && compareAll(Ad.softtype, Bd.softtype, show=show, skip=skip)
  return TP
end

function compareAllSpecial(A::T1,
                           B::T2;
                           skip=Symbol[],
                           show::Bool=true) where {T1 <: GenericFunctionNodeData, T2 <: GenericFunctionNodeData}
  if T1 != T2
    return false
  else
    return compareAll(A, B, skip=skip, show=show)
  end
end

# Compare FunctionNodeData
function compare(a::GenericFunctionNodeData{T1},b::GenericFunctionNodeData{T2}) where {T1, T2}
  # TODO -- beef up this comparison to include the gwp
  TP = true
  TP = TP && a.eliminated == b.eliminated
  TP = TP && a.potentialused == b.potentialused
  TP = TP && a.edgeIDs == b.edgeIDs
  # TP = TP && typeof(a.fnc) == typeof(b.fnc)
  return TP
end

"""
    $SIGNATURES

Compare that all fields are the same in a `::FactorGraph` factor.
"""
function compareFactor(A::DFGFactor,
                       B::DFGFactor;
                       show::Bool=true,
                       skip::Vector{Symbol}=Symbol[],
                       skipsamples::Bool=true,
                       skipcompute::Bool=true  )
  #
  TP =  compareAll(A, B, skip=union([:attributes;:solverData;:_variableOrderSymbols],skip), show=show)
  # TP = TP & compareAll(A.attributes, B.attributes, skip=[:data;], show=show)
  TP = TP & compareAllSpecial(getSolverData(A), getSolverData(B), skip=union([:fnc], skip), show=show)
  TP = TP & compareAllSpecial(getSolverData(A).fnc, getSolverData(B).fnc, skip=union([:cpt;:measurement;:params;:varidx;:threadmodel], skip), show=show)
  TP = TP & (skipsamples || compareAll(getSolverData(A).fnc.measurement, getSolverData(B).fnc.measurement, show=show, skip=skip))
  TP = TP & (skipcompute || compareAll(getSolverData(A).fnc.params, getSolverData(B).fnc.params, show=show, skip=skip))
  TP = TP & (skipcompute || compareAll(getSolverData(A).fnc.varidx, getSolverData(B).fnc.varidx, show=show, skip=skip))

  return TP
end
  # Ad = getSolverData(A)
  # Bd = getSolverData(B)
  # TP =  compareAll(A, B, skip=[:attributes;:data], show=show)
  # TP &= compareAll(A.attributes, B.attributes, skip=[:data;], show=show)
  # TP &= compareAllSpecial(getSolverData(A).fnc, getSolverData(B).fnc, skip=[:cpt;], show=show)
  # TP &= compareAll(getSolverData(A).fnc.cpt, getSolverData(B).fnc.cpt, show=show)

"""
    $SIGNATURES

Compare all variables in both `::FactorGraph`s A and B.

Notes
- A and B should all the same variables and factors.

Related:

`compareFactorGraphs`, `compareSimilarVariables`, `compareVariable`, `ls`
"""
function compareAllVariables(fgA::G1,
                             fgB::G2;
                             skip::Vector{Symbol}=Symbol[],
                             show::Bool=true,
                             skipsamples::Bool=true )::Bool where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  # get all the variables in A or B
  xlA =  listVariables(fgA)
  xlB =  listVariables(fgB)
  vars = union(xlA, xlB)

  # compare all variables exist in both A and B
  TP = length(xlA) == length(xlB)
  for xla in xlA
    TP &= xla in xlB
  end
  # slightly redundant, but repeating opposite direction anyway
  for xlb in xlB
    TP &= xlb in xlA
  end

  # compare each variable is the same in both A and B
  for var in vars
    TP = TP && compareVariable(getVariable(fgA, var), getVariable(fgB, var), skipsamples=skipsamples, skip=skip)
  end

  # return comparison result
  return TP
end

"""
    $SIGNATURES

Compare similar labels between `::FactorGraph`s A and B.

Notes
- At least one variable label should exist in both A and B.

Related:

`compareFactorGraphs`, `compareAllVariables`, `compareSimilarFactors`, `compareVariable`, `ls`.
"""
function compareSimilarVariables(fgA::G1,
                                 fgB::G2;
                                 skip::Vector{Symbol}=Symbol[],
                                 show::Bool=true,
                                 skipsamples::Bool=true )::Bool where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  xlA = listVariables(fgA)
  xlB = listVariables(fgB)

  # find common variables
  xlAB = intersect(xlA, xlB)
  TP = length(xlAB) > 0

  # compare the common set
  for var in xlAB
    @info var
    TP &= compareVariable(getVariable(fgA, var), getVariable(fgB, var), skipsamples=skipsamples, skip=skip)
  end

  # return comparison result
  return TP
end

"""
    $SIGNATURES

Compare similar factors between `::FactorGraph`s A and B.

Related:

`compareFactorGraphs`, `compareSimilarVariables`, `compareAllVariables`, `ls`.
"""
function compareSimilarFactors(fgA::G1,
                               fgB::G2;
                               skipsamples::Bool=true,
                               skipcompute::Bool=true,
                               show::Bool=true  )::Bool where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  xlA = listFactors(fgA)
  xlB = listFactors(fgB)

  # find common variables
  xlAB = intersect(xlA, xlB)
  TP = length(xlAB) > 0

  # compare the common set
  for var in xlAB
    TP = TP && compareFactor(getFactor(fgA, var), getFactor(fgB, var), skipsamples=skipsamples, skipcompute=skipcompute, show=show)
  end

  # return comparison result
  return TP
end

"""
    $SIGNATURES

Compare and return if two factor graph objects are the same, by comparing similar variables and factors.

Notes:
- Default items to skip with `skipsamples`, `skipcompute`.
- User defined fields to skip can be specified with `skip::Vector{Symbol}`.
- To enable debug messages for viewing which fields are not the same:
  - https://stackoverflow.com/questions/53548681/how-to-enable-debugging-messages-in-juno-julia-editor

Related:

`compareSimilarVariables`, `compareSimilarFactors`, `compareAllVariables`, `ls`.
"""
function compareFactorGraphs(fgA::G1,
                             fgB::G2;
                             skipsamples::Bool=true,
                             skipcompute::Bool=true,
                             skip::Vector{Symbol}=Symbol[],
                             show::Bool=true  )::Bool where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  skiplist = Symbol[:g;:bn;:IDs;:fIDs;:id;:nodeIDs;:factorIDs;:fifo;:solverParams; :factorOperationalMemoryType]
  skiplist = union(skiplist, skip)
  @warn "compareFactorGraphs will skip comparisons on: $skiplist"

  TP = compareAll(fgA, fgB, skip=skiplist, show=show)
  TP = TP && compareSimilarVariables(fgA, fgB, skipsamples=skipsamples, show=show, skip=skiplist )
  TP = TP && compareSimilarFactors(fgA, fgB, skipsamples=skipsamples, skipcompute=skipcompute, show=show )
  TP = TP && compareAll(fgA.solverParams, fgB.solverParams, skip=skiplist)

  return TP
end


text:
# push!(LOAD_PATH, "../src")

# using CASIM

# render_option = Not_Render()
# sim_horizon = 5000
# num_run = 3
# init_num_ac = 100
# spawn_rates = [250., 500., 1000., 1500., 2000., 4000., 6000., 8000., 10000.] ./ 3600
# exp = "Simple_10km_init_ac_" * string(init_num_ac)

# """
# 	Simple + Nearest
# """
# cas = Simple_CAS()
# cas_ctrl = Uniform_CAS(cas)
# sensor = Nearest_Intruder_Tracker()
# sensor_ctrl = Uniform_Sensor(sensor)

# for run in 1:num_run
# 	for spawn_rate in spawn_rates
# 		scheme = "Simple_sr_" *
# 			string(Int(spawn_rate * 3600)) * "_run_" * string(run)
# 		airspace_ctrl = Constant_Spawn_Rate_Controller(init_num_ac, spawn_rate)
# 		simulate(exp, scheme, sim_horizon, airspace_ctrl, cas_ctrl, sensor_ctrl, render_option)
# 	end
# end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
#=
testODE:
- Julia version:
- Author: ymocquar
- Date: 2019-11-13
=#
include("../src/interface.jl")
include("../src/henon_heiles.jl")
using LinearAlgebra
using Plots
using Plots.PlotMeasures
function fctmain(n_tau, prec)
    setprecision(prec)
    u0 = BigFloat.([90, -44, 83, 13] // 100)
    tab_eps = zeros(BigFloat, 8)
    epsilon = big"0.19"
    for i = 1:size(tab_eps, 1)
        tab_eps[i] = epsilon
        epsilon /= 10
    end
    nbmaxtest = 12
    order = 6
    t_max = big"1.0"
    y = ones(Float64, nbmaxtest, size(tab_eps, 1))
    x = zeros(Float64, nbmaxtest)
    ind = 1
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]
    for epsilon in tab_eps
        prob = HiOscODEProblem(henon_heiles, u0, (big"0.0", t_max), missing, A, epsilon)
        tabsol = Array{Array{BigFloat,1},1}(undef, 0)
        eps_v = convert(Float32, epsilon)
        nb = 60 * 2^nbmaxtest
        res = solve(prob, nb_tau = n_tau, order = order, nb_t = nb, dense = false)
        par_u0 = res.par_u0
        solRef = res[end]
        nb = 100
        indc = 1
        labels = Array{String,2}(undef, 1, ind)
        while indc <= nbmaxtest
            res = solve(
                prob,
                nb_tau = n_tau,
                order = order,
                nb_t = nb,
                par_u0 = par_u0,
                dense = false,
            )
            sol = res[end]
            nm = norm(sol - solRef, Inf)
            y[indc, ind] = (nm < 1) ? nm : NaN
            x[indc] = 1.0 / nb
            nb *= 2
            indc += 1
        end
        for i = 1:ind
            labels[1, i] = "  ε = $(convert(Float32,tab_eps[i])) "
        end
        p = Plots.plot(
            x,
            view(y, :, 1:ind),
            xlabel = "Δt",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :bottomright,
            label = labels,
            marker = 2,
        )
        pp = Plots.plot(
            x,
            view(y, :, 1:ind),
            xlabel = "Δt",
            xaxis = :log,
            ylabel = "error",
            yaxis = :log,
            legend = :bottomright,
            label = labels,
            marker = 2,
            bottom_margin = 30px,
            left_margin = 60px,
        )
        prec_v = precision(BigFloat)
        Plots.savefig(p, "out/r2p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon_hh_2.pdf")
        Plots.savefig(p, "out/r2p_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon_hh_2.png")
        Plots.savefig(pp, "out/r2pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon_hh_2.pdf")
        Plots.savefig(pp, "out/r2pp_$(prec_v)_$(eps_v)_$(order)_$(n_tau)_epsilon_hh_2.png")
        ind += 1
    end
end

fctmain(32, 256)


text:
using MarketData

facts("constructors") do

  context("inner constructor enforces constraints") do
      @fact 1 => roughly(1.0000000001)
  end

  context("empty constructor") do
      @fact 1 => roughly(1.0000000001)
  end

  context("add entry to existing Blotter object") do
    @fact 1.0000000000000001 => roughly(1)
  end

  context("constructed from orderbook") do
    @fact 1.0000000000000001 => roughly(1)
  end

  context("constructed from signal and FinancialTimeSeries") do
    @fact 1.0000000000000001 => roughly(1)
  end

  context("fillblotter") do
    @fact 1.0000000000000001 => roughly(1)
  end
end

facts("base imports") do

  context("length is correct") do
    @fact_throws log(-1)
    @fact_throws sum(foo, bar)
  end

  context("getindex from single Int") do
    @fact_throws log(-1)
  end

  context("getindex from Int array") do
    @fact_throws log(-1)
  end

  context("getindex from colname") do
    @fact_throws log(-1)
  end

  context("getindex from single date") do
    @fact_throws log(-1)
  end

  context("getindex from date array") do
    @fact_throws log(-1)
  end

  context("getindex from date range") do
    @fact_throws log(-1)
  end
end


text:
module PropertyDestructuringTest

using PropertyDestructuring
using Test

struct Point2D{T}
    x::T
    y::T
end

struct Point3D{T}
    x::T
    y::T
    z::T
end

@testset "PropertyDestructuring.jl" begin

    @testset "Assignment" begin
        pt1 = Point2D(1, 2)
        @destructure (; x, y) = pt1
        @test x == 1
        @test y == 2

        pt2 = Point2D(3.14, 9.99)
        @destructure (; x, y) = pt2
        @test x == 3.14
        @test y == 9.99

        pt3 = Point3D(1.2, 3.4, 5.6)
        @destructure (; x, y) = pt3
        @test x == 1.2
        @test y == 3.4

        @destructure (; x, y) = (x=:x, y="歪")
        @test x === :x
        @test y == "歪"
    end

    @testset "Assignment With Type Annotation" begin
        let pt1 = Point2D(1, 2)
            @destructure (; x::Float64, y::Float64) = pt1
            @test typeof(x) === typeof(y) === Float64
            @test x == 1.0
            @test y == 2.0
        end

        let pt2 = Point2D(3.14, 9.99)
            @destructure (; x::Float64, y::Float64) = pt2
            @test typeof(x) === typeof(y) === Float64
            @test x == 3.14
            @test y == 9.99
        end

        let pt3 = Point3D(1.2f0, 3.4f0, 5.6f0)
            @destructure (; x::Float64, y::Float64) = pt3
            @test typeof(x) === typeof(y) === Float64
            @test x ≈ 1.2 rtol=√eps(Float32)
            @test y ≈ 3.4 rtol=√eps(Float32)
        end

        let pt3 = Point3D(1.2f0, 3.4f0, 5.6f0)
            @destructure (; x::Float32, y::Float32) = pt3
            @test typeof(x) === typeof(y) === Float32
            @test x == 1.2f0
            @test y == 3.4f0
        end
    end

    @testset "Complex Assignment" begin
        pt1 = Point2D(1, 2)
        @destructure (_DUMMY, (; x, y)) = (:DUMMY, pt1)
        @test x == 1
        @test y == 2

        pt2 = Point2D(3.14, 9.99)
        @destructure ((; x, y), _DUMMY) = (pt2, :DUMMY)
        @test x == 3.14
        @test y == 9.99

        pt3 = Point3D(1.2, 3.4, 5.6)
        @destructure ((; x, y), ((; z, w), _DUMMY)) = (pt3, ((z=:z, w="W"), :DUMMY))
        @test x == 1.2
        @test y == 3.4
        @test z === :z
        @test w == "W"

        let pt3 = Point3D(1.2, 3.4, 5.6), other = (a=1.1, b=2.2)
            @destructure ((; x, y), (; a::Float16, b::Float32)) = (pt3, other)
            typeof(a) === Float16
            typeof(b) === Float32
            @test x == 1.2
            @test y == 3.4
            @test a ≈ Float16(1.1)
            @test b ≈ 2.2f0
        end
    end

    @testset "Function Argument (1)" begin
        @destructure function _test1((; x, y))
            (x, y)
        end

        pt1 = Point2D(1, 2)
        x, y = _test1(pt1)
        @test x == 1
        @test y == 2

        pt2 = Point2D(3.14, 9.99)
        x, y = _test1(pt2)
        @test x == 3.14
        @test y == 9.99

        pt3 = Point3D(1.2, 3.4, 5.6)
        x, y = _test1(pt3)
        @test x == 1.2
        @test y == 3.4

        x, y = _test1((x=:x, y="歪"))
        @test x === :x
        @test y == "歪"
    end

    @testset "Function Argument (2)" begin
        @destructure _test2((; x, y)) = (x, y)

        pt1 = Point2D(1, 2)
        x, y = _test2(pt1)
        @test x == 1
        @test y == 2

        pt2 = Point2D(3.14, 9.99)
        x, y = _test2(pt2)
        @test x == 3.14
        @test y == 9.99

        pt3 = Point3D(1.2, 3.4, 5.6)
        x, y = _test2(pt3)
        @test x == 1.2
        @test y == 3.4

        x, y = _test2((x=:x, y="歪"))
        @test x === :x
        @test y == "歪"

        @destructure myreal((; re)::Complex) = re
        c = rand(ComplexF64)
        @test myreal(c) == real(c) == c.re
    end

    @testset "`in` clause of `for`" begin
        pts = [
            Point2D(1, 2),
            Point2D(3.14, 9.99),
            Point3D(1.2, 3.4, 5.6),
            (x=:x, y="歪"),
        ]

        index = 1
        @destructure for (; x, y) in pts
            @test x == pts[index].x
            @test y == pts[index].y
            index += 1
        end

        index = 1
        @destructure for (; x::Float64, y::Float64) in pts[1:3]
            @test typeof(x) === typeof(y) === Float64
            @test x == pts[index].x
            @test y == pts[index].y
            index += 1
        end

        @destructure for (pt, (; x, y)) in zip(pts, pts)
            @test x == pt.x
            @test y == pt.y
        end

        index = 1
        @destructure for (; x, y) in pts, _DUMMY in 1
            @test x == pts[index].x
            @test y == pts[index].y
            index += 1
        end

        @destructure for (pt, (; x, y)) in zip(pts, pts), _DUMMY in 1
            @test x == pt.x
            @test y == pt.y
        end

        @destructure for (pt, (; x::Float64, y::Float64)) in zip(pts[1:3], pts)
            @test typeof(x) === typeof(y) === Float64
            @test x == pt.x
            @test y == pt.y
        end
    end

    @testset "Combined" begin
        let
            @destructure check_combined((;x, y), x_expected, y_expected) = x == x_expected && y == y_expected

            pts = [
                Point2D(1, 2),
                Point2D(3.14, 9.99),
                Point3D(1.2, 3.4, 5.6),
                (x=:x, y="歪"),
            ]

            @destructure for (pt, (; x, y)) in zip(pts, pts)
                @test check_combined(pt, x, y)
            end

            for pt in pts
                @destructure (;x, y) = pt
                @test check_combined(pt, x, y)
            end
        end

        @destructure let
            check_combined((;x, y), x_expected, y_expected) = x == x_expected && y == y_expected

            pts = [
                Point2D(1, 2),
                Point2D(3.14, 9.99),
                Point3D(1.2, 3.4, 5.6),
                (x=:x, y="歪"),
            ]

            for (pt, (; x, y)) in zip(pts, pts)
                @test check_combined(pt, x, y)
            end

            for pt in pts
                (;x, y) = pt
                @test check_combined(pt, x, y)
            end
        end
    end

    @testset "Const" begin
        include("consttest.jl")

        @test typeof(Test20211215.num) === UInt8
        @test typeof(Test20211215.den) === UInt8
        @test Test20211215.num == 0x01
        @test Test20211215.den == 0x05
    end
end

end  # module


text:
#Test functional Integrate and adjoint for euler equation.

@doc """
Euler Equation -- test_adjoint

The function tests for the correctness of objective and non-objective functional
computation and then tests if the adjoint vector is being computed correctly.
Adjoint vector computation is based on an objective function. This is a serial
test and uses the input file called

`input_vals_vortex_adjoint_DG.jl`

"""->

function test_adjoint()

  @testset "Checking complete derivative of a functional using adjoint vector" begin
    # println("testing adjoint functions\n")
    #fname = "input_vals_vortex_adjoint_DG.jl"
    fname = "input_vals_airfoil.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    @assert opts["aoa"] == 2.0

    lift = createFunctional(mesh, sbp, eqn, opts, 1)

    h = 1e-20
    pert = Complex128(0, h)

    pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts)
    ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
    adjoint_vec = zeros(Complex128, mesh.numDof)
    calcAdjoint(mesh, sbp, eqn, opts, ls, lift, adjoint_vec, recalc_jac=true, recalc_pc=true)

    # compute dR/daoa and dJ/daoa
    eqn.params.aoa += pert
    evalResidual(mesh, sbp, eqn, opts)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
    dRdaoa = imag(eqn.res_vec)/h

    val = evalFunctional(mesh, sbp, eqn, opts, lift)
    dJdaoa = imag(val)/h
    eqn.params.aoa -= pert

    DJDaoa = dJdaoa + dot(adjoint_vec, dRdaoa)

    # save value for parallel tests
    fname = "BoundaryForceData_values.dat"
    file = open(fname, "w")
    println(file, val)
    close(file)

    # Check complete derivatives w.r.t alpha using finite difference
    pert = 1e-6
    eqn.params.aoa += pert
    solvePDE(mesh, sbp, eqn, opts, mesh)
    val2 = evalFunctional(mesh, sbp, eqn, opts, lift)
    eqn.params.aoa -= pert

    DJDaoa_fd = (val2 - real(val))/pert
    err_val = norm(DJDaoa - DJDaoa_fd, 2)

    println("dJDaoa = ", DJDaoa)
    println("dJdaoa_fd = ", DJDaoa_fd)
    println("diff = ", DJDaoa - DJDaoa_fd)
    @test isapprox( err_val, 0.0) atol= 1e-6

  end # End  testset("Checking complete derivative of a functional using adjoint vector")

end # End function test_adjoint

add_func1!(EulerTests, test_adjoint, [TAG_ADJOINT, TAG_LONGTEST])


text:
"""
    timespan(eph::Ephem)

This function returns the first and last time available in the ephemeris file associated to eph.

# Arguments:
- `eph` : ephemeris

# Return:
a tuple containing:
    * firsttime: Julian date of the first time
    * lasttime: Julian date of the last time
    * continuous: information about the availability of the quantities over the time span

        It returns the following value in the parameter continuous :

        1 if the quantities of all bodies are available for any time between the first and last time.
        2 if the quantities of some bodies are available on discontinuous time intervals between the first and last time.
        3 if the quantities of each body are available on a continuous time interval between the first and last time,
          but not available for any time between the first and last time.

See: https://www.imcce.fr/content/medias/recherche/equipes/asd/calceph/html/c/calceph.multiple.html#menu-calceph-gettimespan
"""
function timespan(eph::Ephem)
    @_checkPointer eph.data "Ephemeris is not properly initialized!"

    firsttime = Ref{Cdouble}(0)
    lasttime  = Ref{Cdouble}(0)
    continous = Ref{Cint}(0)
    stat = ccall((:calceph_gettimespan, libcalceph), Cint, (Ptr{Cvoid},Ref{Cdouble},Ref{Cdouble},Ref{Cint}),
                                                            eph.data,firsttime,lasttime,continous)

    @_checkStatus stat "Unable to compute ephemeris"
    return firsttime[],lasttime[],continous[]
end


text:
# --------------------------------------------------------------------------- #
# Loading an instance of SPP (format: OR-library)

function loadSPP(fname)
    f=open(fname)
    # lecture du nbre de contraintes (m) et de variables (n)
    m, n = parse.(Int, split(readline(f)) )
    # lecture des n coefficients de la fonction economique et cree le vecteur d'entiers c
    C = parse.(Int, split(readline(f)) )
    # lecture des m contraintes et reconstruction de la matrice binaire A
    A=zeros(Int, m, n)
    for i=1:m
        # lecture du nombre d'elements non nuls sur la contrainte i (non utilise)
        readline(f)
        # lecture des indices des elements non nuls sur la contrainte i
        for valeur in split(readline(f))
          j = parse(Int, valeur)
          A[i,j]=1
        end
    end
    close(f)
    return C, A
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 4fc69b90-72c2-11eb-3c8d-5ffa667de01f
using GMT, CSV, DataFrames, HTTP

# ╔═╡ 91c16f80-72c1-11eb-0d19-5903f6ef7d1e
md"
# Choropleth maps
A choropleth is a thematic map where areas are colored in proportion to a variable such as population density. GMT lets us plot choropleth maps but the process is not straightforward because we need to put the color information in the header of a multi-segment file. To facilitate this, some tools were added to the Julia wrapper.
"

# ╔═╡ 55e740a0-72c3-11eb-3680-370bee35effa
md"
Load packages needed to download data and put into a DataFrame
"

# ╔═╡ 5b77eca0-72c2-11eb-161c-6b6c3581cede
paises = CSV.File(HTTP.get("https://raw.githubusercontent.com/tillnagel/unfolding/master/data/data/countries-population-density.csv").body, delim=';') |> DataFrame;

# ╔═╡ eb759ed0-72c5-11eb-16b6-191c48be50da
md"
Extract the european countries from the DCW file
"

# ╔═╡ 63e05440-72c2-11eb-1c2f-851a0cf3b553
# Extract the Europe countries from the DCW file
eu = coast(DCW="=EU+z", dump=true);

# ╔═╡ f4fb64b0-72c2-11eb-18d2-9720e67e3dad
md"
Generate two vector with country codes (2 chars codes) and the population density for European countries
"

# ╔═╡ 82ec2cb0-72c2-11eb-1140-95c37ebdef78
codes, vals = GMT.mk_codes_values(paises[!, 2], paises[!, 3], region="eu");

# ╔═╡ a5fa0560-72c2-11eb-3b9b-950358050b24
md"
Create a Categorical CPT that allows plotting the choropleth.
Note that we need to limit the CPT range to leave out the very high density *states* like Monaco otherwise all others would get the same color.
"

# ╔═╡ b9626f6e-72c2-11eb-0bd3-b3fc48195d9e
C = cpt4dcw(codes, vals, range=[10, 500, 10]);

# ╔═╡ 02f51f20-72c3-11eb-02e7-436644e34088
md"
Show the result
"

# ╔═╡ 1e4c05e2-72c3-11eb-32a2-89200f2e26d1
plot(eu, cmap=C, fill="+z", proj=:guess, region="-76/36/26/84",
	title="Population density", colorbar=true, show=true)

# ╔═╡ Cell order:
# ╟─91c16f80-72c1-11eb-0d19-5903f6ef7d1e
# ╟─55e740a0-72c3-11eb-3680-370bee35effa
# ╠═4fc69b90-72c2-11eb-3c8d-5ffa667de01f
# ╠═5b77eca0-72c2-11eb-161c-6b6c3581cede
# ╟─eb759ed0-72c5-11eb-16b6-191c48be50da
# ╠═63e05440-72c2-11eb-1c2f-851a0cf3b553
# ╟─f4fb64b0-72c2-11eb-18d2-9720e67e3dad
# ╠═82ec2cb0-72c2-11eb-1140-95c37ebdef78
# ╟─a5fa0560-72c2-11eb-3b9b-950358050b24
# ╠═b9626f6e-72c2-11eb-0bd3-b3fc48195d9e
# ╟─02f51f20-72c3-11eb-02e7-436644e34088
# ╠═1e4c05e2-72c3-11eb-32a2-89200f2e26d1


text:
module DL

using Printf
using LinearAlgebra
import Base.*, Base.+

export (*), (+), ReLU, VariableTensor, MSE, ConstantTensor, backward, zero_grad, SGD, Linear, Tensor

# Tensor
abstract type Tensor end
abstract type LeafTensor <: Tensor end
abstract type FunctionalTensor <: Tensor end

#= 掛け算 t = t1 * t2 =#
mutable struct MulTensor <: FunctionalTensor
    x
    tensors::Array{Tensor, 1}
end
(*)(t3::Tensor, t4::Tensor) = abstract_operator(MulTensor, (x1, x2) -> x1 * x2, [t3, t4])
# t2 <- t1 <- t3 * t4
function grad_fn(t1::MulTensor, grad_up)
    t3, t4 = t1.tensors
    # dt2/dt1 * dt1/dt3 = grad_up * t4.x
    # dt2/dt1 * dt1/dt4 = grad_up * t3.x
    return [t4.x * grad_up, grad_up * t3.x]
end

#= 足し算 =#
mutable struct AddTensor <: FunctionalTensor
    x
    tensors::Array{Tensor, 1}
end
# forward
(+)(t3::Tensor, t4::Tensor) = abstract_operator(AddTensor, (x1, x2) -> x1 .+ x2, [t3, t4])
# backward t2 <- t1 <- t3 + t4
function grad_fn(t1::AddTensor, grad_up)
    t3, t4 = t1.tensors
    # dt2/dt1 * dt1/dt3 = grad_up * I
    # dt2/dt1 * dt1/dt4 = grad_up * I
    return [grad_up, grad_up]
end

#= ReLU function t = ReLU(t1) =#
mutable struct ReLUTensor <: FunctionalTensor
    x
    tensors::Array{Tensor, 1}
end
# forward
function ReLU(t::Tensor)
    return abstract_operator(ReLUTensor, x -> x .* (x .> 0), [t])
end
# backward t2 <- t1 <- ReLU(t3)
function grad_fn(t1::ReLUTensor, grad_up)
    t3 = t1.tensors[1]
    # dt2/dt1 * dt1/dt3 = grad_up * d ReLU(x)/dx
    return [grad_up .* (t3.x' .> 0)]
end

#= Mean Squared Error t = MSE(t1, t2) =#
mutable struct MSETensor <: FunctionalTensor
    x
    tensors::Array{Tensor, 1}
end
# forward
function MSE(t1::Tensor, t2::Tensor)
    mse = (x1, x2) -> sum((x1 - x2).^2) / length(x1)
    return abstract_operator(MSETensor, mse, [t1, t2])
end
# t2 <- t1 <- MSE(t3, t4)
function grad_fn(t1::MSETensor, grad_up)
    t3, t4 = t1.tensors
    # dt2/dt1 * dt1/dt3 = grad_up * d/dx sum((x - y)^2) / length(x)
    N = length(t3.x)
    dt1dt3 = 2*(t3.x - t4.x)' / N
    dt1dt4 = 2*(t4.x - t3.x)' / N
    return [grad_up * dt1dt3, grad_up * dt1dt4]
end

#= t2 <- t1 末端の部分 微分あり =#
mutable struct VariableTensor <: LeafTensor
    x
    grad
    VariableTensor(x) = new(x, zero(x))
end

#= ただの定数 =#
mutable struct ConstantTensor <: LeafTensor
    x
end

# backward 計算
# 最初
function backward(t::Tensor)
	backward(t, 1)
end
# 再帰
function backward(t::Tensor, grad_up)
    # 微分はtが持つTensorsの数だけあるから，リストで返すことにする．
    vgrad = grad_fn(t, grad_up)
	# 微分と対応するTensorのペアをつくる
    pair = zip(t.tensors, vgrad)
    # 掘り下げ
	map(x -> backward(x[1], x[2]), pair)
end
# ConstantTensorとか
function backward(t::LeafTensor, grad_up)
end
# VariableTensorは勾配を蓄積していく
function backward(t::VariableTensor, grad_up)
    t.grad += grad_up'
end

# 抽象化された関数
# 演算子を計算するときにConstant同士の演算の出力はConstantであるようにする
# Array{T, 1} where T <: Tensor という書き方でTensor以下のタイプのTのArrayを網羅できる
function abstract_operator(EnzansiTensor::DataType, f::Function, tensors::Array{T, 1} where T <: Tensor)
    # Tensor.xを展開してfに代入する
    y = f(map(tensor -> tensor.x, tensors)...)
    return EnzansiTensor(y, tensors)
end
function abstract_operator(EnzansiTensor::DataType, f::Function, tensors::Array{ConstantTensor, 1})
    y = f(map(tensor -> tensor.x, tensors)...)
    return ConstantTensor(y)
end

function tensor_map(f::Function, t::Tensor)
    #= 任意の関数をVariableTensorに作用させる
    毎回掘り下げる必要はないと思うので考えなおす余地あり
    =#
    # 掘り下げていく
    map(t -> tensor_map(f, t), t.tensors)
end
function tensor_map(f::Function, t::LeafTensor)
    # 末端だけどVariableTensor以外のtypeはパラメータを変えない
    nothing
end
function tensor_map(f::Function, t::VariableTensor)
    # VariableTensorにだけ作用させる
    f(t)
end

function zero_grad(tensor::Tensor)
    f = t -> t.grad = zero(t.grad)
    tensor_map(f, tensor)
end

function SGD(lr=1e-3)
    f = t::Tensor -> t.x = t.x - lr * t.grad
    optimizer = t::Tensor -> tensor_map(f, t)
    return optimizer
end

function Linear(in_features::Int64, out_features::Int64)
    W = VariableTensor(randn(out_features, in_features))
    b = VariableTensor(zeros(out_features, 1))
    @show b.x
    return x -> W*x + b
end

function printarray(x)
    show(stdout, "text/plain", x)
    @printf "\n"
end

end

max_stars_repo_name: Soyukke/julia-ayame-machine-learning

max_issues_repo_name: Soyukke/julia-ayame-machine-learning

max_forks_repo_name: Soyukke/julia-ayame-machine-learning


text:
using Flux: Dropout, rand!, _dropout_kernel

_dropout_shape(s, dims...) = tuple((i ∈ dims ? 1 : si for (i, si) ∈ enumerate(s))...)

function (a::Dropout)(x, dims)
  a.active || return x
  y = similar(x, _dropout_shape(size(x), dims...)...)
  rand!(y)
  y .= _dropout_kernel.(y, a.p, 1 - a.p)
  return x .* y
end