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

using Markdown
using InteractiveUtils

# ╔═╡ 90bdce42-7980-11eb-08ad-0f39e65e1864
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

# ╔═╡ 20f35fc2-799a-11eb-12e9-4588c897199d
include("utils.jl");

# ╔═╡ 51dfdce2-7980-11eb-229a-3b3196ded703
md"""
## C03/w03: Identifying safe loans with decision trees

The [LendingClub](https://www.lendingclub.com/) is a peer-to-peer leading company that directly connects borrowers and potential lenders/investors. In this notebook, we will build a classification model to predict whether or not a loan provided by LendingClub is likely to [default](https://en.wikipedia.org/wiki/Default_%28finance%29).

In this notebook we will use data from the LendingClub to predict whether a loan will be paid off in full or the loan will be [charged off](https://en.wikipedia.org/wiki/Charge-off) and possibly go into default. In this assignment we will:

  - Use SFrames to do some feature engineering.
  - Train a decision-tree on the LendingClub dataset.
  - Predict whether a loan will default along with prediction probabilities (on a validation set).
  - Train a complex tree model and compare it to simple tree model.
"""

# ╔═╡ a36cad10-7980-11eb-212f-67774be5fd35
md"""
### Load LendingClub dataset

We will be using a dataset from the [LendingClub](https://www.lendingclub.com/). A parsed and cleaned form of the dataset is availiable [here](https://github.com/learnml/machine-learning-specialization-private). Make sure you **download the dataset** before running the following command.
"""

# ╔═╡ a36b25ee-7980-11eb-2308-e1b576d9883e
begin
  loans =  CSV.File(("../../ML_UW_Spec/C03/data/lending-club-data.csv");
  header=true) |> DataFrame;
  first(loans, 3)
end

# ╔═╡ a36b0642-7980-11eb-3a6f-e3f48bc58616
md"""
### Exploring some features

Let's quickly explore what the dataset looks like. First, let's print out the column names to see what features we have in this dataset.

"""

# ╔═╡ 19b49852-7981-11eb-1d65-99a8f9119422
length(names(loans)), names(loans)

# ╔═╡ 882c144c-7989-11eb-0e5d-4d2b0d13b459
function by_feature(df::DataFrame, feature::Symbol)
  groupby(df, feature) |>
    df_ -> combine(df_, nrow => Symbol(feature, "_count"))
end

# ╔═╡ e84e8f1c-798e-11eb-23e6-0dafc1deaff4
function barplot(df::DataFrame, feature::Symbol)
  let df=by_feature(df, feature), fcnt=Symbol(feature, "_count"), yaxe=df[!, feature]

    bar(df[!, fcnt], yticks=(1:length(yaxe), collect(yaxe)),
      orientation=:h, label=String(feature))
  end
end

# ╔═╡ d5331906-798f-11eb-21c5-bf17a158c639
barplot(loans, :grade)

# ╔═╡ 87ede65a-798d-11eb-2610-7db61052073b
md"""
We can see that over half of the loan grades are assigned values B or C. Each loan is assigned one of these grades, along with a more finely discretized feature called sub_grade (feel free to explore that feature column as well!). These values depend on the loan application and credit report, and determine the interest rate of the loan. More information can be found here.

Now, let's look at a different feature.
"""

# ╔═╡ 9405f590-798d-11eb-3a55-619e7d2c5cb6
barplot(loans, :sub_grade)

# ╔═╡ 93e75248-798d-11eb-2e4b-cb8072c75df7
barplot(loans, :home_ownership)

# ╔═╡ 93ce629a-798d-11eb-02fe-7b0556848660
md"""
This feature describes whether the loanee is mortaging, renting, or owns a home. We can see that a small percentage of the loanees own a home.
"""

# ╔═╡ 93b3c4aa-798d-11eb-04e9-f526070ac3eb
md"""
### Exploring the target column

The target column (label column) of the dataset that we are interested in is called `bad_loans`. In this column **1** means a *risky* (bad) loan **0** means a *safe* loan.

In order to make this more intuitive and consistent with the lectures, we reassign the target to be:
  - **+1** as a safe  loan,
  - **-1** as a risky (bad) loan.

We put this in a new column called `safe_loans`.

"""

# ╔═╡ 939cde84-798d-11eb-244a-9ba95ee1ada7
insertcols!(loans, :safe_loans => ifelse.(loans.bad_loans .== 0, 1, -1),
  makeunique=true);

# ╔═╡ 937f41c6-798d-11eb-0eec-ffd7cb58b2b6
md"""
Now, let us explore the distribution of the column safe_loans. This gives us a sense of how many safe and risky loans are present in the dataset.
"""

# ╔═╡ 9368145e-798d-11eb-2232-dd7572dfb755
barplot(loans, :safe_loans)

# ╔═╡ 9350b5f4-798d-11eb-2c9f-a1cbd7526d1a
with_terminal() do
  num_loans = size(loans)[1]

  println((num_safe_loans=round(sum(loans.safe_loans .== 1) / num_loans; digits=4), num_bad_loans=round(sum(loans.safe_loans .== -1) / num_loans; digits=4)))
end

# ╔═╡ 93194734-798d-11eb-3d96-1f8a46085885
md"""
We should have:
  - Around 81% safe loans
  - Around 19% risky loans

It looks like most of these loans are safe loans (thankfully). But this does make our problem of identifying risky loans challenging.
"""

# ╔═╡ 0120efaa-7992-11eb-313e-6b5e4dfd4800
md"""
### Features for the classification algorithm

In this assignment, we will be using a subset of features (categorical and numeric). The features we will be using are described in the code comments below. If you are a finance geek, the LendingClub website has a lot more details about these features.

"""

# ╔═╡ 1671f0ac-7992-11eb-20ff-912bf4353a24
begin
  const Features = [:grade,       # grade of the loan
            :sub_grade,             # sub-grade of the loan
            :short_emp,             # one year or less of employment
            :emp_length_num,        # number of years of employment
            :home_ownership,        # home_ownership status: own, mortgage or rent
            :dti,                   # debt to income ratio
            :purpose,               # the purpose of the loan
            :term,                  # the term of the loan
            :last_delinq_none,      # has borrower had a delinquincy
            :last_major_derog_none, # has borrower had 90 day or worse rating
            :revol_util,            # percent of available credit being used
            :total_rec_late_fee,    # total late fees received to day
           ]

  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

# ╔═╡ 16564a14-7992-11eb-3600-cd4cf1573933
md"""
What remains now is a subset of features and the target that we will use for the rest of this notebook.
"""

# ╔═╡ 0aa86afe-79ff-11eb-0d39-4fee0ee8c7fc
md"""
#### MLJ Data Preparation

MLJ distinguished between (refer to [Data interpretation: Scientific Types](https://alan-turing-institute.github.io/DataScienceTutorials.jl/data/scitype/):
 - encoding or *machine type* (like an Integer, a String....) and
 - how data should be interpreted or *scientific type*

Let's have a look at the current mapping machine type/scientific type on our dataframe.
"""

# ╔═╡ 34c2bb96-79aa-11eb-1d64-d7e7ff02159c
schema(loans)

# ╔═╡ 67b8b802-79fa-11eb-241a-17e808148f7c
## which attribute/feature is Textual?

# map(attr -> (attr, elscitype(loans[!, attr])), Features) |>
#   ary -> filter(((_a, st) = t) -> st == Textual, ary)

# ╔═╡ e9c1bcc8-79fe-11eb-1437-9530a1eefbdb
md"""
Those textual attributes need to be coerced to more meaningful scitype.

Based on the values of those attributes we can infer the folloing mapping:
"""

# ╔═╡ 1d6247c6-7a01-11eb-39e1-c1565af8a588
begin
  struct F_ST
    names::Vector{Symbol}
    st_map::Dict{Symbol, DataType}

    function F_ST(names::Vector{Symbol}, st_map::Dict{Symbol, DataType})
      @assert length(names) ≥ 1
      @assert length(st_map) ≥ 0
      new(names, st_map)
    end
  end

  ## subset of features whose scientific types need to be coerced
  const Features_ST = Dict{Symbol, DataType}(
    :grade => OrderedFactor{7},          # because 7 ≠ values
    :sub_grade => OrderedFactor{35},
    :home_ownership => Multiclass{4},
    :purpose => Multiclass{12},
    :term => OrderedFactor{2},
  )

  features = F_ST(Features, Features_ST)

  const Target_ST = Dict{Symbol, DataType}(
    Target => Multiclass{2},
  )

  target = F_ST([Target], Target_ST)
end

# ╔═╡ f2c7be06-7992-11eb-2c95-979c0965fe77
md"""
### Sample data to balance classes

As we explored above, our data is disproportionally full of safe loans. Let's create two datasets: one with just the safe loans (safe_loans_raw) and one with just the risky loans (risky_loans_raw).

"""

# ╔═╡ f29fef50-7992-11eb-2ed8-2bf706c9828a
begin
  safe_loans_raw = loans[loans[!, target.names[1]] .== 1, :]
  risky_loans_raw = loans[loans[!, target.names[1]] .== -1, :]

  with_terminal() do
    @printf("Number of safe loans  : %d\n", size(safe_loans_raw)[1])
    @printf("Number of risky loans : %d\n", size(risky_loans_raw)[1])
  end
end

# ╔═╡ 1629b576-7992-11eb-1fef-a3538166d17d
md"""
Now, write some code to compute below the percentage of safe and risky loans in the dataset and validate these numbers against what was given earlier in the assignment:
"""

# ╔═╡ 16055156-7992-11eb-3461-17b098dc3a27
begin
  function loan_ratios(safe_loans, risky_loans; n_loans=size(loans)[1])
      (size(safe_loans)[1] / n_loans, size(risky_loans)[1] / n_loans)
  end

  p_safe_loans, p_risky_loans = loan_ratios(safe_loans_raw, risky_loans_raw)

  with_terminal() do
    @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

# ╔═╡ 00f5c508-7992-11eb-1087-e7c63d454e06
md"""
One way to combat class imbalance is to undersample the larger class until the class distribution is approximately half and half. Here, we will undersample the larger class (safe loans) in order to balance out our dataset. This means we are throwing away many data points.

We used seed=42.
"""

# ╔═╡ 00d20af2-7992-11eb-188e-55d5bf3fc1b3
begin
  ## Since there are fewer risky loans than safe loans, find the ratio of the sizes
  ## and use that percentage to undersample the safe loans.
  nsl = size(safe_loans_raw)[1]
  perc = size(risky_loans_raw)[1] / nsl

  risky_loans = risky_loans_raw
  ## down sample safe loans
  safe_loans = safe_loans_raw[sampling(nsl, perc; seed=42), :]
  loans_data = [risky_loans; safe_loans]
end

# ╔═╡ 1f4430b2-7999-11eb-1fd1-ddf50fb7a6dd
md"""
Now, let's verify that the resulting percentage of safe and risky loans are each nearly 50%.
"""

# ╔═╡ 1f29904a-7999-11eb-1007-4da27f7c6677
begin
  np_safe_loans, np_risky_loans = loan_ratios(safe_loans, risky_loans, n_loans=size(loans_data)[1])

  with_terminal() do
    @printf("Percentage of safe loans: %1.2f%%\n", np_safe_loans * 100.)
    @printf("Percentage of risky loans: %1.2f%%\n", np_risky_loans * 100.)
    @printf("Total number of loans in our new dataset: %d\n", size(loans_data)[1])
  end
end

# ╔═╡ 1f0ac002-7999-11eb-0bb5-63faa4dcb6cc
md"""
**Note**: There are many approaches for dealing with imbalanced data, including some where we modify the learning algorithm. These approaches are beyond the scope of this course.

For this assignment, we use the simplest possible approach, where we subsample the overly represented class to get a more balanced dataset. In general, and especially when the data is highly imbalanced, we recommend using more advanced methods.
"""

# ╔═╡ f0a6c596-7999-11eb-0629-979335f89922
md"""
### Split data into training and validation sets

We split the data into training and validation sets using an 80/20 split and specifying `seed=1` so everyone gets the same results.

**Note**: In previous assignments, we have called this a **train-test split**. However, the portion of data that we don't train on will be used to help **select model parameters** (this is known as model selection). Thus, this portion of data should be called a **validation set**. Recall that examining performance of various potential models (i.e. models with different parameters) should be on validation set, while evaluation of the final selected model should always be on test data.

Typically, we would also save a portion of the data (a real test set) to test our final model on or use cross-validation on the training set to select our final model. But for the learning purposes of this assignment, we won't do that.
"""

# ╔═╡ f0890d6e-7999-11eb-3482-4769d6489475
begin
  train_data, validation_data = train_test_split(loans_data; split=0.8, seed=42);
  size(train_data), size(validation_data)
end

# ╔═╡ 68e89cf0-79fb-11eb-24ff-d37b75bebe2e
first(train_data, 3)

# ╔═╡ f070b8f2-7999-11eb-2849-f56a3c74659b
md"""
### Use decision tree to build a classifier

Now, let's use the MLJ decision tree learner to create a loan prediction model on the training data. (In the next assignment, you will implement your own decision tree learning algorithm.)

Our feature columns and target column have already been decided above.

"""

# ╔═╡ a7f44f2a-799c-11eb-34ff-9b3613ca036c
# Pkg.add("MLJDecisionTreeInterface")

# ╔═╡ 971c2ac2-799c-11eb-1d64-d7e7ff02159c
@load DecisionTreeClassifier pkg=DecisionTree

# ╔═╡ 77751a08-79a0-11eb-2c16-f342efd58086
function dt_runner(df::DF;
    features=features, target=target,
    max_depth=-1) where {DF <: Union{DataFrame, SubDataFrame}}
  X, y = data_prep_Xy(df, features, target)

  # dt_clf = MLJDecisionTreeInterface.DecisionTreeClassifier(
  #   max_depth = max_depth,
  #   min_samples_leaf = 1,
  #   min_samples_split = 2,
  #   min_purity_increase = 0.0,
  #   n_subfeatures = 0,
  #   post_prune = false,
  #   merge_purity_threshold = 1.0,
  #   pdf_smoothing = 0.0,
  #   display_depth = 5
  # )

  pipeclf = @pipeline(
    OneHotEncoder(),
    MLJDecisionTreeInterface.DecisionTreeClassifier()
  )

  pipeclf.decision_tree_classifier.max_depth=max_depth

  dt_mach = machine(pipeclf, X, y)
  fit!(dt_mach, verbosity=0)
  dt_mach
end

# ╔═╡ 92fe3f40-798d-11eb-190c-81209cbb835b
# function dt_runner_alt(df::DF;
#     features=features, target=traget
#     max_depth=3) where {DF <: Union{DataFrame, SubDataFrame}}
#   X, y = data_prep_Xy(df, features, target)

#   pipeclf = @pipeline(
#     OneHotEncoder(),
#     MLJDecisionTreeInterface.DecisionTreeClassifier()
#   )

#   r_mpi = range(pipeclf, :(pipeclf.decision_tree_classifier.max_depth),
#     lower=1, upper=max_depth)
#   # r_msl = range(pipeclf, :(pipeclf.decision_tree_classifier.min_samples_leaf),
#   #   lower=1, upper=50)

#   tm = TunedModel(model=pipeclf, ranges=[r_mpi], tuning=Grid(resolution=8),
#     resampling=CV(nfolds=5, rng=112),
#     operation=predict_mode,
#     measure=misclassification_rate,
#     check_measure=false)

#   mtm = machine(tm, Xt, y)
#   fit!(mtm, rows=1:size(df)[1], verbosity=1)

#   mtm
# end

# ╔═╡ 08fc6234-799e-11eb-072e-f31cc84a982a
begin
  dt_model = dt_runner(train_data; max_depth=6)

  with_terminal() do
    println("Summary dt_model:")
    println("--------------------------------------")
    println(fitted_params(dt_model))
    println("--------------------------------------")
  end
end

# ╔═╡ f76d8430-79a7-11eb-1d64-d7e7ff02159c
md"""

Error / References:

  - https://github.com/alan-turing-institute/MLJ.jl/issues/423

  - https://stackoverflow.com/questions/66365311/the-scitype-of-x-in-machine-is-incompatible-with-model

  - https://alan-turing-institute.github.io/MLJ.jl/stable/common_mlj_workflows/#Data-ingestion-1

"""

# ╔═╡ dcc73932-79ac-11eb-196a-4d0b50beb96d
with_terminal() do
  println(report(dt_model))
end

# ╔═╡ d8d0090a-79a0-11eb-1d64-d7e7ff02159c
md"""
#### Building a smaller tree

A tree can be hard to visualize graphically, and moreover, it may overfit.. Here, we instead learn a smaller model with max depth of 2 to gain some intuition and to understand the learned tree more.
"""

# ╔═╡ 6f5b75f0-7a03-11eb-1e21-ab19398446c2
begin
  dt_small_model = dt_runner(train_data; max_depth=2)

  with_terminal() do
    println("Summary dt_small_model:")
    println("--------------------------------------")
    println(fitted_params(dt_small_model))
    println("--------------------------------------")
  end
end

# ╔═╡ 0bd3a464-7a18-11eb-2e7d-c7246d8cc30d
md"""
#### Building a bigger tree
"""

# ╔═╡ 1f101802-7a18-11eb-3e07-e5f88b6f9b58
begin
  dt_big_model = dt_runner(train_data; max_depth=10)

  with_terminal() do
    println("Summary dt_big_model:")
    println("--------------------------------------")
    println(fitted_params(dt_big_model))
    println("--------------------------------------")
  end
end

# ╔═╡ 6f3bc75a-7a03-11eb-1df3-392d8e3e0363
md"""
#### Making predictions

Let's consider two positive and two negative examples from the validation set and see what the model predicts. We will do the following:

  - Predict whether or not a loan is safe.
  - Predict the probability that a loan is safe.

"""

# ╔═╡ 6f2014a6-7a03-11eb-1ec1-abd58d12aac7
begin
  n_validation_data = data_prep(validation_data, features)

  valid_safe_loans = filter(target.names[1] => ==(1), n_validation_data)
  valid_risky_loans = filter(target.names[1] => ==(-1), n_validation_data)

  sample_valid_data = vcat(valid_safe_loans[1:2, :], valid_risky_loans[1:2, :])
  sample_valid_data;
end

# ╔═╡ c6f5a1ea-7a0b-11eb-0dd9-f79fc57abf0e
begin
  using Distributions

  y = sample_valid_data.safe_loans
  ŷ = predict(dt_model, select(sample_valid_data, Features))

  L = levels(y)
  broadcast(pdf, ŷ, y), pdf(ŷ, L)
end

# ╔═╡ a0f6600c-7a08-11eb-0190-a1eea0f3ef62
function get_levels(col)
  length(levels(col)), levels(col)
end

# ╔═╡ 3ff06776-7a08-11eb-0187-4191b4637cbe
# length(levels(train_data.grade)), levels(train_data.grade)
get_levels(train_data.grade)

# ╔═╡ 747f7110-7a08-11eb-16d1-6333aa1c0df5
get_levels(train_data.sub_grade)

# ╔═╡ 74601010-7a08-11eb-164b-8de382485329
get_levels(sample_valid_data.grade)

# ╔═╡ b39d761c-7a05-11eb-2f0f-27f32329cb2f
md"""
#### Explore label predictions

Now, we will use our model  to predict whether or not a loan is likely to default. For each row in the `sample_validation_data`, use the `dt_model` to predict whether or not the loan is classified as a **safe loan**.
"""

# ╔═╡ b3833b4e-7a05-11eb-09f6-43a205481b3b
begin
  ŷ_mode = predict_mode(dt_model,
    select(sample_valid_data, Features))
  println(ŷ_mode)

  perc_correct_preds = sum(ŷ_mode .== sample_valid_data.safe_loans) / size(sample_valid_data)[1]
  (percentage_prediction=perc_correct_preds,)
end

# ╔═╡ b368f158-7a05-11eb-17e9-8ddaea0b2386
md"""
**Quiz Question:** What percentage of the predictions on `sample_valid_data` did `dt_model` get correct?
  - cf. above cell
"""

# ╔═╡ d927b594-7a0c-11eb-3329-11481535ab85
begin
  ix = argmax(broadcast(pdf, ŷ, y))
  ŷ_prob = broadcast(pdf, ŷ, y)
  (ix=ix, value=ŷ_prob[ix])
end

# ╔═╡ b34f5fb8-7a05-11eb-1399-d9cd7c8ded35
md"""
**Quiz Question**: Which loan has the highest probability of being classified as a safe loan?
  - in this case the third (index 3, *Julia indexes start from 1*)

"""

# ╔═╡ 1194ab16-7a15-11eb-127e-1dec3928fc2a
md"""
**Checkpoint:** Can you verify that for all the predictions with `probability >= 0.5`, the model predicted the label **+1**?
"""

# ╔═╡ b8c7a3c0-7a1d-11eb-296e-ef4d312ffcfa
ŷ_prob, y

# ╔═╡ c89eeb2e-7a1c-11eb-2278-ab37447cf22d
## Checkpoint code

# begin
#   ## indexes for which y is 1
#   ixes_ = filter(((_, v)=t) -> v,
#       collect(enumerate(sample_valid_data.safe_loans .== 1))) |>
#     a_ -> map(t -> t[1], a_)

#   ## all probs matching y .== 1 are ≥ 0.5
#   @test all(ŷ_prob[ixes_, :] .≥ 0.5)
# end

# ╔═╡ 47718414-7a1f-11eb-0b31-6d187e952bad
md"""
##### Tricky predictions!

Now, we will explore something pretty interesting. For each row in the `sample_valid_data`, what is the probability (according to dt_small_model) of a loan being classified as safe?

"""

# ╔═╡ c850fd10-7a1c-11eb-3b93-1bbeb7098830
begin
  # ŷ_vec = predict_mode(dt_small_model, select(sample_valid_data, Features))
  ŷ_dt_small_ = predict(dt_small_model, select(sample_valid_data, Features))

  # ŷ_vec,
  broadcast(pdf, ŷ_dt_small_, y)
end

# ╔═╡ 12909320-7a15-11eb-13df-691a76b3a823
md"""
### Evaluating accuracy of the decision tree model

Recall that the accuracy is defined as follows:

```math
\mbox{accuracy} = \frac{\mbox{\# correctly classified examples}}{\mbox{\# total examples}}
```

Let us start by evaluating the accuracy of the `dt_small_model` and `dt_model` on the training data
"""

# ╔═╡ 931f5666-7a15-11eb-0042-cd16bc34e020
function accuracy(ŷ, y)
  sum(ŷ .== y) / size(y)[1]
end

# ╔═╡ b3352738-7a05-11eb-13eb-f341e0e2e6df
begin
  ys = n_validation_data.safe_loans

  ŷ_dt = predict_mode(dt_model, select(n_validation_data, Features))
  ŷ_dt_small = predict_mode(dt_small_model, select(n_validation_data, Features))
  ŷ_dt_big = predict_mode(dt_big_model, select(n_validation_data, Features))

  # mode.(y_hat)
  # broadcast(pdf, y_hats, ys, ) #[1:2]
  acc_dt = accuracy(ŷ_dt, ys)
  acc_dt_small = accuracy(ŷ_dt_small, ys)
  acc_dt_big = accuracy(ŷ_dt_big, ys)

  (acc_dt=acc_dt |> a_ -> round(a_; digits=3),
    acc_dt_small=acc_dt_small |> a_ -> round(a_; digits=3),
    acc_dt_big=acc_dt_big |> a_ -> round(a_; digits=3))
end

# ╔═╡ 4dce4096-7a19-11eb-0bd8-f14ec24ed701
md"""
### Quantifying the cost of mistakes

Every mistake the model makes costs money. In this section, we will try and quantify the cost of each mistake made by the model.

Assume the following:

* **False negatives**: Loans that were actually safe but were predicted to be risky. This results in an oppurtunity cost of losing a loan that would have otherwise been accepted.
* **False positives**: Loans that were actually risky but were predicted to be safe. These are much more expensive because it results in a risky loan being given.
* **Correct predictions**: All correct predictions don't typically incur any cost.

Let's write code that can compute the cost of mistakes made by the model. Complete the following 4 steps:
1. First, let us compute the predictions made by the model.
1. Second, compute the number of false positives.
2. Third, compute the number of false negatives.
3. Finally, compute the cost of mistakes made by the model by adding up the costs of true positives and false positives.

"""

# ╔═╡ fdd913ca-7a1a-11eb-3bee-6f04974e3141
function calc_fp(ŷ, y)
  filter(((ŷ_, y_)=t) -> ŷ_ == y_ == 1,
    collect(zip(ŷ .== 1, y .== -1))) |> length
end

# ╔═╡ c111df1c-7a1a-11eb-1474-b96dcad21315
function calc_fn(ŷ, y)
  filter(((ŷ_, y_)=t) -> ŷ_ == y_ == 1,
    collect(zip(ŷ .== -1, y .== 1))) |> length
end

# ╔═╡ 66bc0b4c-7a19-11eb-1d94-df9fa1b7e7b5
md"""
False positives are predictions where the model predicts +1 but the true label is -1. Complete the following code block for the number of false positives:
"""

# ╔═╡ 7120d40c-7a1a-11eb-33f2-8f6b79fcb83a
fp_dt, fp_dt_big = calc_fp(ŷ_dt, ys),  calc_fp(ŷ_dt_big, ys)

# ╔═╡ 6d3396cc-7a19-11eb-0c3f-3d0456bff8c5
md"""
False negatives are predictions where the model predicts -1 but the true label is +1. Complete the following code block for the number of false negatives:
"""

# ╔═╡ 6d176dec-7a19-11eb-25ae-09cf985eb22a
fn_dt, fn_dt_big = calc_fn(ŷ_dt, ys), calc_fn(ŷ_dt_big, ys)

# ╔═╡ 6cf9fa46-7a19-11eb-249a-d300a87184cf
md"""
**Quiz Question:** Let us assume that each mistake costs money:
  - Assume a cost of \$10,000 per false negative.
  - Assume a cost of \$20,000 per false positive.

What is the total cost of mistakes made by `decision_tree_model` on `validation_data`?
"""

# ╔═╡ 8b58e726-7a1a-11eb-076b-e3b226df3fc2
(dt_model_cost=fp_dt * 20000 + fn_dt * 10000, dt_big_model_cost=fp_dt_big * 20000 + fn_dt_big * 10000)

# ╔═╡ Cell order:
# ╟─51dfdce2-7980-11eb-229a-3b3196ded703
# ╠═90bdce42-7980-11eb-08ad-0f39e65e1864
# ╟─a36cad10-7980-11eb-212f-67774be5fd35
# ╠═a36b25ee-7980-11eb-2308-e1b576d9883e
# ╟─a36b0642-7980-11eb-3a6f-e3f48bc58616
# ╠═19b49852-7981-11eb-1d65-99a8f9119422
# ╠═882c144c-7989-11eb-0e5d-4d2b0d13b459
# ╠═e84e8f1c-798e-11eb-23e6-0dafc1deaff4
# ╠═d5331906-798f-11eb-21c5-bf17a158c639
# ╟─87ede65a-798d-11eb-2610-7db61052073b
# ╠═9405f590-798d-11eb-3a55-619e7d2c5cb6
# ╠═93e75248-798d-11eb-2e4b-cb8072c75df7
# ╟─93ce629a-798d-11eb-02fe-7b0556848660
# ╟─93b3c4aa-798d-11eb-04e9-f526070ac3eb
# ╠═939cde84-798d-11eb-244a-9ba95ee1ada7
# ╟─937f41c6-798d-11eb-0eec-ffd7cb58b2b6
# ╠═9368145e-798d-11eb-2232-dd7572dfb755
# ╠═9350b5f4-798d-11eb-2c9f-a1cbd7526d1a
# ╟─93194734-798d-11eb-3d96-1f8a46085885
# ╟─0120efaa-7992-11eb-313e-6b5e4dfd4800
# ╠═1671f0ac-7992-11eb-20ff-912bf4353a24
# ╟─16564a14-7992-11eb-3600-cd4cf1573933
# ╟─0aa86afe-79ff-11eb-0d39-4fee0ee8c7fc
# ╠═34c2bb96-79aa-11eb-1d64-d7e7ff02159c
# ╠═67b8b802-79fa-11eb-241a-17e808148f7c
# ╟─e9c1bcc8-79fe-11eb-1437-9530a1eefbdb
# ╠═1d6247c6-7a01-11eb-39e1-c1565af8a588
# ╟─f2c7be06-7992-11eb-2c95-979c0965fe77
# ╠═f29fef50-7992-11eb-2ed8-2bf706c9828a
# ╟─1629b576-7992-11eb-1fef-a3538166d17d
# ╠═16055156-7992-11eb-3461-17b098dc3a27
# ╟─00f5c508-7992-11eb-1087-e7c63d454e06
# ╠═00d20af2-7992-11eb-188e-55d5bf3fc1b3
# ╟─1f4430b2-7999-11eb-1fd1-ddf50fb7a6dd
# ╠═1f29904a-7999-11eb-1007-4da27f7c6677
# ╟─1f0ac002-7999-11eb-0bb5-63faa4dcb6cc
# ╟─f0a6c596-7999-11eb-0629-979335f89922
# ╠═20f35fc2-799a-11eb-12e9-4588c897199d
# ╠═f0890d6e-7999-11eb-3482-4769d6489475
# ╠═68e89cf0-79fb-11eb-24ff-d37b75bebe2e
# ╟─f070b8f2-7999-11eb-2849-f56a3c74659b
# ╠═a7f44f2a-799c-11eb-34ff-9b3613ca036c
# ╠═971c2ac2-799c-11eb-1d64-d7e7ff02159c
# ╠═77751a08-79a0-11eb-2c16-f342efd58086
# ╠═92fe3f40-798d-11eb-190c-81209cbb835b
# ╠═08fc6234-799e-11eb-072e-f31cc84a982a
# ╟─f76d8430-79a7-11eb-1d64-d7e7ff02159c
# ╠═dcc73932-79ac-11eb-196a-4d0b50beb96d
# ╟─d8d0090a-79a0-11eb-1d64-d7e7ff02159c
# ╠═6f5b75f0-7a03-11eb-1e21-ab19398446c2
# ╟─0bd3a464-7a18-11eb-2e7d-c7246d8cc30d
# ╠═1f101802-7a18-11eb-3e07-e5f88b6f9b58
# ╟─6f3bc75a-7a03-11eb-1df3-392d8e3e0363
# ╠═6f2014a6-7a03-11eb-1ec1-abd58d12aac7
# ╠═a0f6600c-7a08-11eb-0190-a1eea0f3ef62
# ╠═3ff06776-7a08-11eb-0187-4191b4637cbe
# ╠═747f7110-7a08-11eb-16d1-6333aa1c0df5
# ╠═74601010-7a08-11eb-164b-8de382485329
# ╟─b39d761c-7a05-11eb-2f0f-27f32329cb2f
# ╠═b3833b4e-7a05-11eb-09f6-43a205481b3b
# ╟─b368f158-7a05-11eb-17e9-8ddaea0b2386
# ╠═c6f5a1ea-7a0b-11eb-0dd9-f79fc57abf0e
# ╠═d927b594-7a0c-11eb-3329-11481535ab85
# ╟─b34f5fb8-7a05-11eb-1399-d9cd7c8ded35
# ╟─1194ab16-7a15-11eb-127e-1dec3928fc2a
# ╠═b8c7a3c0-7a1d-11eb-296e-ef4d312ffcfa
# ╠═c89eeb2e-7a1c-11eb-2278-ab37447cf22d
# ╟─47718414-7a1f-11eb-0b31-6d187e952bad
# ╠═c850fd10-7a1c-11eb-3b93-1bbeb7098830
# ╟─12909320-7a15-11eb-13df-691a76b3a823
# ╠═931f5666-7a15-11eb-0042-cd16bc34e020
# ╠═b3352738-7a05-11eb-13eb-f341e0e2e6df
# ╟─4dce4096-7a19-11eb-0bd8-f14ec24ed701
# ╠═fdd913ca-7a1a-11eb-3bee-6f04974e3141
# ╠═c111df1c-7a1a-11eb-1474-b96dcad21315
# ╟─66bc0b4c-7a19-11eb-1d94-df9fa1b7e7b5
# ╠═7120d40c-7a1a-11eb-33f2-8f6b79fcb83a
# ╟─6d3396cc-7a19-11eb-0c3f-3d0456bff8c5
# ╠═6d176dec-7a19-11eb-25ae-09cf985eb22a
# ╟─6cf9fa46-7a19-11eb-249a-d300a87184cf
# ╠═8b58e726-7a1a-11eb-076b-e3b226df3fc2


text:
function get_boundary_coordinates_from_index(bndy_index::Int64, qtEl::QuadTreeMeshes.QuadTreeElement, qt::QuadTreeMeshes.QuadTree)
  # get quad tree element

  if bndy_index <= 4
    b1, b2 = qtEl.bbLeftBottomIndex, qtEl.bbRightBottomIndex
  elseif bndy_index <= 8
    b1, b2 = qtEl.bbRightBottomIndex, qtEl.bbRightTopIndex
  elseif bndy_index <= 12
    b1, b2 = qtEl.bbRightTopIndex, qtEl.bbLeftTopIndex
  else
    b1, b2 = qtEl.bbLeftTopIndex, qtEl.bbLeftBottomIndex
  end

  inner_index = (bndy_index - 1) % 4
  x1, x2 = qt.vertices[b1], qt.vertices[b2]
  ds = x2 - x1
  pos = x1 + ds * (0.125 + inner_index * 0.25)

  return pos
end

function check_leave_intersection(qt::QuadTreeMeshes.QuadTree, elIndex::QuadTreeMeshes.ElIndex, ls::GeometryTypes.Simplex{2, QuadTreeMeshes.Point}, center_vertex::Nullable{QuadTreeMeshes.vertex_index})
  qt_element = qt.elements[elIndex]
  mesh_element = get(qt.values[elIndex])

  tinterPoints = Array{QuadTreeMeshes.Point, 1}()

  sb1, sb2 = qt.vertices[qt_element.bbLeftBottomIndex], qt.vertices[qt_element.bbRightBottomIndex]
  sb1_intersects, sb1i = GeometryTypes.intersects(ls, GeometryTypes.LineSegment(sb1, sb2))
  eb1, eb2 = qt.vertices[qt_element.bbRightBottomIndex], qt.vertices[qt_element.bbRightTopIndex]
  eb1_intersects, eb1i = GeometryTypes.intersects(ls, GeometryTypes.LineSegment(eb1, eb2))
  nb1, nb2 = qt.vertices[qt_element.bbRightTopIndex], qt.vertices[qt_element.bbLeftTopIndex]
  nb1_intersects, nb1i = GeometryTypes.intersects(ls, GeometryTypes.LineSegment(nb1, nb2))
  wb1, wb2 = qt.vertices[qt_element.bbLeftBottomIndex], qt.vertices[qt_element.bbLeftTopIndex]
  wb1_intersects, wb1i = GeometryTypes.intersects(ls, GeometryTypes.LineSegment(wb1, wb2))
  if sb1_intersects
    push!(tinterPoints, sb1i)
  end
  if eb1_intersects
    push!(tinterPoints, eb1i)
  end
  if nb1_intersects
    push!(tinterPoints, nb1i)
  end
  if wb1_intersects
    push!(tinterPoints, wb1i)
  end

  # remove duplicates
  interPoints = Array{QuadTreeMeshes.Point, 1}()
  for tp in tinterPoints
    found = false
    for p in interPoints
      if !found && dot((p-tp), (p-tp)) < 1e-10
        found = true
      end
    end
    if !found
      push!(interPoints, tp)
    end
  end

  # sort according to distance from start
  sort!(interPoints, by = p -> dot((p-ls[1]), (p-ls[1])))

  npoints = length(interPoints)

  center_is_in_element = false
  if !isnull(center_vertex)
    vPoint = qt.vertices[get(center_vertex)]
    vpx, vpy = vPoint
    if ((vpx >= sb1[1] && vpx <= sb2[1]) && (vpy >= eb1[2] && vpy <= eb2[2]))
      # vertex is in element
      center_is_in_element = true
      @test(!isnull(mesh_element.center) && get(mesh_element.center) == get(center_vertex))
    else
      @test(isnull(mesh_element.center))
    end
  end
  if center_is_in_element
    # we must have at most two boundaries - one that starts at the center
    # and one that ends there (there might have been one deleted already though)
    @assert(length(mesh_element.boundaries) <= 2, "Element $elIndex is $mesh_element.")
    s, e = ls
    ls_starts_at_vertex = dot(s-vPoint, s-vPoint) <= 1e-10
    ls_ends_at_vertex = dot(e-vPoint, e-vPoint) <= 1e-10
    @assert((ls_starts_at_vertex && !ls_ends_at_vertex) || (!ls_starts_at_vertex && ls_ends_at_vertex))
    found = Nullable{Tuple{Int, QuadTreeMeshes.BoundaryVertex}}()
    for (bidx, b) in enumerate(mesh_element.boundaries)
      if ls_starts_at_vertex && b[1].vt == QuadTreeMeshes.Inner && b[1].vertex == get(center_vertex)
        found = Nullable{Tuple{Int, QuadTreeMeshes.BoundaryVertex}}((bidx, b[2]))
      elseif ls_ends_at_vertex && b[2].vt == QuadTreeMeshes.Inner && b[2].vertex == get(center_vertex)
        found = Nullable{Tuple{Int, QuadTreeMeshes.BoundaryVertex}}((bidx, b[1]))
      end
    end
    @assert(!isnull(found))
    foundBoundary, foundVertex = get(found)
    boundary_vertex_point = qt.vertices[foundVertex.vertex]
    @assert(length(interPoints) == 1)
    dist = boundary_vertex_point - interPoints[1]
    @test(dot(dist, dist) <= 1e-10)
    mesh_element.boundaries = deleteat!(mesh_element.boundaries, foundBoundary)
  else
    @assert(length(interPoints) <= 2)
    if length(interPoints) > 0
      if (length(interPoints) == 1)
        input_boundary = (interPoints[1], interPoints[1])
      elseif length(interPoints) == 2
        input_boundary = (interPoints[1], interPoints[2])
      end
      found = Nullable{Int}()
      for (ind, b) in enumerate(mesh_element.boundaries)
        b1, b2 = b
        b1p = qt.vertices[b1.vertex]
        b2p = qt.vertices[b2.vertex]
        db1 = b1p - input_boundary[1]
        db2 = b2p - input_boundary[2]
        if dot(db1, db1) <= 1e-10 && dot(db2, db2) <= 1e-10
          found = Nullable{Int}(ind)
        end
      end
      if isnull(found)
        print_with_color(:red, "No suitable boundaries found. Intersection points $(interPoints)\n")
        for b in mesh_element.boundaries
          b1, b2 = b
          print_with_color(:green, "$(b1.vertex) [$(qt.vertices[b1.vertex])] -> $(b2.vertex) [$(qt.vertices[b2.vertex])]\n")
        end
      end
      @assert(!isnull(found))
      mesh_element.boundaries = deleteat!(mesh_element.boundaries, get(found))
    end
  end
end

function subdividefunc(x)
  #println("$x")
end

function onlinesegment(
  xi::QuadTreeMeshes.Point,
  x1::QuadTreeMeshes.Point,
  x2::QuadTreeMeshes.Point)
    if (xi ≈ x1) || (xi ≈ x2)
      return true
    end
    local d1i = normalize(xi - x1)
    local di2 = normalize(x2 - xi)
    return d1i ≈ di2
end

function orientation(v1::QuadTreeMeshes.Point, v2::QuadTreeMeshes.Point, v3::QuadTreeMeshes.Point)
  local xa = v1[1]
  local ya = v1[2]
  local xb = v2[1]
  local yb = v2[2]
  local xc = v3[1]
  local yc = v3[2]

  local det = (xa*yb - xa*yc - xb*ya + xb*yc + xc*ya - xc*yb)
  if abs(det) <= 1e-10
    return 0
  elseif det > 0.
    return 1
  else
    return 2
  end
end
"""
    pointsintersect(p1::QuadTreeMeshes.Point, q1::QuadTreeMeshes.Point, p2::QuadTreeMeshes.Point, q2::QuadTreeMeshes.Point)

Checks if line segments given by their end points intersect.

# Remarks
* Taken from geeksforgeeks.org/check-if-two-given-line-segments-intersect
"""
function pointsintersect(p1::QuadTreeMeshes.Point, q1::QuadTreeMeshes.Point, p2::QuadTreeMeshes.Point, q2::QuadTreeMeshes.Point)
  local o1 = orientation(p1,q1,p2)
  local o2 = orientation(p1,q1,q2)
  local o3 = orientation(p2,q2,p1)
  local o4 = orientation(p2,q2,q1)
  if (o1 != o2 && o3 != o4)
    return true
  elseif (o1 == 0 && onlinesegment(p2, p1, q1))
    return true
  elseif (o2 == 0 && onlinesegment(q2, p1, q1))
    return true
  elseif (o3 == 0 && onlinesegment(p1, p2, q2))
    return true
  elseif (o4 == 0 && onlinesegment(q1, p2, q2))
    return true
  else
    return false
  end
end

function line_intersects_line(ls1::GeometryTypes.LineSegment{QuadTreeMeshes.Point}, ls2::GeometryTypes.LineSegment{QuadTreeMeshes.Point})
  return pointsintersect(ls1[1], ls1[2], ls2[1], ls2[2])
end

function line_intersects_rectangle(r::GeometryTypes.SimpleRectangle{Float64}, ls::GeometryTypes.LineSegment{QuadTreeMeshes.Point})
  ((x1, y1), (x2, y2)) = ls
  xbl, ybl, xtr, ytr = r.x, r.y, r.x + r.w, r.y + r.h
  f(x, y) = (y2 - y1) * x + (x1 - x2) * y + (x2 * y1 - x1 * y2)
  s1 = f(xbl, ybl)
  s2 = f(xbl, ytr)
  s3 = f(xtr, ybl)
  s4 = f(xtr, ytr)
  if (s1>0 && s2>0 && s3>0 && s4>0) || (s1<0 && s2<0 && s3<0 && s4<0)
    return false
  end
  if (x1 > xtr && x2 > xtr)
    return false
  end
  if (x1 < xbl && x2 < xbl)
    return false
  end
  if (y1 > ytr && y2 > ytr)
    return false
  end
  if (y1 < ybl && y2 < ybl)
    return false
  end
  return true
end

function subdivide_random!(qt::QuadTreeMeshes.QuadTree)
  elIndex = rand(1:size(qt.elements, 1))
  while (QuadTreeMeshes.has_child(qt, elIndex))
    elIndex = rand(1:size(qt.elements, 1))
  end

  QuadTreeMeshes.subdivide!(qt, elIndex, subdividefunc)
end

function get_neighbour_from_index(qt::QuadTreeMeshes.QuadTree, elIndex::Int, neighIndex::Int)
  el = qt.elements[elIndex]
  if (neighIndex == 1)
    return el.northWest
  end
  if (neighIndex == 2)
    return el.northEast
  end
  if (neighIndex == 3)
    return el.southWest
  end
  if (neighIndex == 4)
    return el.southEast
  end
end

function check_subdivision_levels(qt::QuadTreeMeshes.QuadTree)
  segments = Array{GeometryTypes.LineSegment{GeometryTypes.Point{2, Float64}},1}()

  # collect all boundary edges of leaves
  for curNode in qt.elements
    if (isnull(curNode.northWest))
      bb = QuadTreeMeshes.get_element_bounding_box(qt, curNode)
      push!(segments, GeometryTypes.LineSegment(GeometryTypes.Point{2, Float64}(bb.x, bb.y), GeometryTypes.Point{2,Float64}(bb.x + bb.w, bb.y)))
      push!(segments, GeometryTypes.LineSegment(GeometryTypes.Point{2, Float64}(bb.x + bb.w, bb.y), GeometryTypes.Point{2, Float64}(bb.x + bb.w, bb.y + bb.h)))
      push!(segments, GeometryTypes.LineSegment(GeometryTypes.Point{2, Float64}(bb.x, bb.y), GeometryTypes.Point{2, Float64}(bb.x, bb.y + bb.h)))
      push!(segments, GeometryTypes.LineSegment(GeometryTypes.Point{2, Float64}(bb.x, bb.y + bb.h), GeometryTypes.Point{2, Float64}(bb.x + bb.w, bb.y + bb.h)))
    end
  end

  # intersecting boundary edges need to have at least one endpoint in common
  while !isempty(segments)
    curSeg = pop!(segments)

    for testSeg in segments
      r0, r1 = testSeg
      t0, t1 = curSeg
      # check if they are parallel and instersect
      if (dot((r0 - r1), (t0-t1))) != 0
        intersects = line_intersects_line(curSeg, testSeg)
        if intersects
          @test isapprox(r0, t0) || isapprox(r0, t1) || isapprox(r1, t0) || isapprox(r1, t1)
          if !(isapprox(r0, t0) || isapprox(r0, t1) || isapprox(r1, t0) || isapprox(r1, t1))
            @assert(false)
          end
        end
      end
    end
  end
end

ext: jl


text:
# This file is a part of SimilaritySearch.jl

export BinaryHammingDistance

"""
   BinaryHammingDistance()

Binary hamming uses bit wise operations to count the differences between bit strings
"""
struct BinaryHammingDistance <: PreMetric end

"""
    evaluate(::BinaryHammingDistance, a, b)
    evaluate(::BinaryHammingDistance, a::AbstractVector, b::AbstractVector) where {T<:Unsigned}

Computes the binary hamming distance for bit types and arrays of bit types
"""
function evaluate(::BinaryHammingDistance, a, b)
    d = 0
    @inbounds @simd for i in eachindex(a)
        d += count_ones(a[i] ⊻ b[i])
    end

    d
end

function evaluate(::BinaryHammingDistance, a::T, b::T)::Float64 where {T<:Unsigned}
    count_ones(a ⊻ b)
end


text:
using Push
using Base.Test

cfg_path = joinpath(dirname(@__FILE__), "configuration/name.cfg")
cfg = Push.load_configuration(cfg_path)

# NAME.QUOTE
s = Push.run("(NAME.QUOTE NAME.QUOTE)", cfg)
@test s.name == [convert(Symbol, "NAME.QUOTE")]

# NAME.=
s = Push.run("(X Y NAME.=)", cfg)
@test s.boolean == [false]
s = Push.run("(X X NAME.=)", cfg)
@test s.boolean == [true]

# NAME.DUP
s = Push.run("(X NAME.DUP)", cfg)
@test s.name == [:X, :X]

# NAME.FLUSH
s = Push.run("(X Y Z NAME.FLUSH)", cfg)
@test isempty(s.name)

# NAME.POP
s = Push.run("(X Y Z NAME.POP)", cfg)
@test s.name == [:X, :Y]

# NAME.ROT
s = Push.run("(Y Z NAME.ROT)", cfg)
@test s.name == [:Y, :Z]
s = Push.run("(A X Y Z NAME.ROT)", cfg)
@test s.name == [:A, :Z, :X, :Y]

# NAME.SHOVE
s = Push.run("(A B C D 0 NAME.SHOVE)", cfg)
@test s.name == [:A, :B, :C, :D]
s = Push.run("(A B C D -100 NAME.SHOVE)", cfg)
@test s.name == [:A, :B, :C, :D]
s = Push.run("(A B C D 1 NAME.SHOVE)", cfg)
@test s.name == [:A, :B, :D, :C]
s = Push.run("(A B C D 2 NAME.SHOVE)", cfg)
@test s.name == [:A, :D, :B, :C]
s = Push.run("(A B C D 3 NAME.SHOVE)", cfg)
@test s.name == [:D, :A, :B, :C]
s = Push.run("(A B C D 10 NAME.SHOVE)", cfg)
@test s.name == [:D, :A, :B, :C]

# NAME.STACKDEPTH
s = Push.run("(A B C D NAME.STACKDEPTH)", cfg)
@test s.name == [:A, :B, :C, :D] && s.integer == [4]
s = Push.run("(NAME.STACKDEPTH)", cfg)
@test s.integer == [0]

# NAME.SWAP
s = Push.run("(A B C NAME.SWAP)", cfg)
@test s.name == [:A, :C, :B]
s = Push.run("(A NAME.SWAP)", cfg)
@test s.name == [:A]

# NAME.YANK
s = Push.run("(A B C D 0 NAME.YANK)", cfg)
@test s.name == [:A, :B, :C, :D]
s = Push.run("(A B C D -100 NAME.YANK)", cfg)
@test s.name == [:A, :B, :C, :D]
s = Push.run("(A B C D 1 NAME.YANK)", cfg)
@test s.name == [:A, :B, :D, :C]
s = Push.run("(A B C D 2 NAME.YANK)", cfg)
@test s.name == [:A, :C, :D, :B]
s = Push.run("(A B C D 3 NAME.YANK)", cfg)
@test s.name == [:B, :C, :D, :A]
s = Push.run("(A B C D 10 NAME.YANK)", cfg)
@test s.name == [:B, :C, :D, :A]

# NAME.YANKDUP
s = Push.run("(A B C D 0 NAME.YANKDUP)", cfg)
@test s.name == [:A, :B, :C, :D, :D]
s = Push.run("(A B C D -100 NAME.YANKDUP)", cfg)
@test s.name == [:A, :B, :C, :D, :D]
s = Push.run("(A B C D 1 NAME.YANKDUP)", cfg)
@test s.name == [:A, :B, :C, :D, :C]
s = Push.run("(A B C D 2 NAME.YANKDUP)", cfg)
@test s.name == [:A, :B, :C, :D, :B]
s = Push.run("(A B C D 3 NAME.YANKDUP)", cfg)
@test s.name == [:A, :B, :C, :D, :A]
s = Push.run("(A B C D 10 NAME.YANKDUP)", cfg)
@test s.name == [:A, :B, :C, :D, :A]

# NAME.RAND
s = Push.run("(NAME.RAND)", cfg)
@test length(s.name) == 1

# NAME.RANDBOUNDNAME
s = Push.run("(NAME.RANDBOUNDNAME)", cfg)
@test length(s.name) == 1 && in(s.name[1], Push.list_instructions(s))


text:
using Documenter, Hecke, Pkg

include(normpath(joinpath(Hecke.pkgdir, "docs", "Build.jl")))

Build.make(Hecke; strict=false, local_build=false, doctest=true, format = :mkdocs)

## Overwrite printing to make the header not full of redundant nonsense
## Turns
##   Hecke.Order - Method
## into
##   Order - Method
#
## To remove the '-'
## Documenter.Utilities.print_signature(io::IO, signature)        = print(io, signature)
#
## To remove the "Method", "Type", "Module" use the following
## Documenter.Utilities.doccat(b::Base.Docs.Binding, ::Type)  = ""
## doccat(::Type)     = ""
## doccat(::Module)   = ""
#
## Remove the module prefix
#Base.print(io::IO, b::Base.Docs.Binding) = print(io, b.var)

# We use 'mike', a mkdocs extension to display a selector.
# We create the versions.json for the selector Usually Documenter.jl creates a
# versions.js containing the information.
function Documenter.Writers.HTMLWriter.generate_version_file(versionfile::AbstractString, entries, symlinks = [])
    versionfile = joinpath(dirname(versionfile), "versions.json")
    open(versionfile, "w") do buf
        symlinks = Dict(symlinks)
        isempty(entries) && return
        println(buf, "[")
        for (j, folder) in enumerate(entries)
          print(buf, "{\"version\": \"$(folder)\", \"title\": \"$(folder)\", \"aliases\": ")
          print(buf, "[")
          for (i, al) in enumerate([e for e in symlinks if e.first == folder])
            print(buf, "\"", al.second, "\"")
            if i < length([e for e in symlinks if e.first == folder])
              print(buf, ", ")
            end
          end
          print(buf, "]}")
          if j < length(entries)
            print(buf, ",")
          end
          print(buf, "\n")
        end
        print(buf, "]")
    end
end

deps = Pkg.dependencies()
ver = Pkg.dependencies()[Base.UUID("3e1990a7-5d81-5526-99ce-9ba3ff248f21")]
s = string(ver.version)
if ENV["RELEASE_VERSION"] == "master"
  s = "dev"
  cmd = `mike deploy $s`
else
  s = ENV["RELEASE_VERSION"]
  cmd = `mike deploy $s`
end

deploydocs(
  repo = "github.com/thofma/Hecke.jl.git",
  deps = Deps.pip("pymdown-extensions",
                  "pygments",
                  "mkdocs",
                  "python-markdown-math",
                  "mkdocs-material",
                  "mkdocs-cinder",
                  "mike"),
  target = "site",
  make = () -> run(cmd),
)


text:
"""
A key table stores the data as a set of dictionaries for each field.
"""
immutable KeyTable{Index,StorageTypes,Key,KeyType} <: AbstractTable{Index,Key}
    data::Dict{KeyType,StorageTypes}

    function KeyTable(data_in)
        check_table(Index,StorageTypes)
        check_key(Key,KeyType)
        new(data_in)
    end
end

@generated check_key{Key,KeyType}(::Key,::Type{KeyType}) = (eltype(Key()) == KeyType) ? (return nothing) : (str = "KeyType $KeyType doesn't match Key $Key"; return :(error(str)))


text:
function get_median_step_infected(df_slice::DataFrame)
    cols = [:pc_infected_avg, :pc_infected_q05, :pc_infected_q95, :pc_infected_lo, :pc_infected_hi]

    res = Pair{Symbol,Int64}[]
    for col in cols
        v = minimum( i -> (abs(0.5 - df_slice[i, col]), df_slice.step[i]), 1:nrow(df_slice))[2]
        push!(res, Symbol(string("step_",col)) => v )
    end
    push!(res, :left_value => 0.0)
    push!(res, :max_step=>maximum(df_slice.step))
    DataFrame(  res...)
end

sigm(x) = exp(x)/(exp(x)+1)

vals = sigm.(-8:0.1:8) #161
med = findfirst(>=(0.5), vals)
1:med

sum(vals[1:med])
sum(1-vals[i] for i in med:length(vals))


text:
# Unit Test for `SolveQuadraticProgram()`

using LinearAlgebra;
using SparseArrays;
using Random;
using StableRNGs;
using Dates;
using BenchmarkTools;
using Printf;
using CSV;
using Tables;
# using DataFrames;

seedNumber = 1234;
Random.seed!(seedNumber);
# Random.seed!(StableRNG(seedNumber), seedNumber);
Random.default_rng() = StableRNG(seedNumber);

include("GenerateQuadraticProgram.jl");
include("SolveQuadraticProgram.jl");
include("LinearSystemSolvers.jl");

# @enum ProblemClass randomQp = 1 inequalityConstrainedQp equalityConstrainedQp optimalControl portfolioOptimization lassoOptimization huberFitting supportVectorMachine isotonicRegression

## Parameters

# Simulaion
numSimulations  = 5;
mNumElements = [    0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;
0010 0100;] .÷ 5;

mNumConstraints = [ 0000 0000;
0000 0000;
0000 0000;
0000 0000;
0000 0000;
0000 0000;
0000 0000;
0000 0000;
0000 0000;]; #<! Set to 0 for OSQP Paper dimensions

benchmarkCsvFileName = "QPSBenchmark.csv";

# Solver
numIterations   = 50000; #<! Enough for Primal Dual / ADMM Convergence
numItrPolish    = 0;
hLinSolInit     = FacLdlInit;
hLinSol!        = FacLdl!;

hSolverFun(vX, mP, vQ, mA, vL, vU) = SolveQuadraticProgram!(vX, mP, vQ, mA, vL, vU, hLinSolInit, hLinSol!; numIterations = numIterations, numItrPolish = numItrPolish);

solverLabel = "QPS LDLFactorizations.jl Solver";
codeVersion = "0.1.000";
cpuInfo     = rstrip(lstrip(Sys.cpu_info()[1].model));
benchDate   = Dates.format(now(Dates.UTC), "yyyy_mm_dd_SS_MM_HH");

# Benchmarking
benchMarkSamples    = 15;
benchMarkEvals      = 1;
benchMarkSeconds    = 120;

numTests    = numSimulations * length(mNumElements);
vTestData   = Vector{Any}(undef, 4 * numTests + 4);
vTestHeader = Vector{Any}(undef, 4 * numTests + 4);

vTestData[1] = solverLabel;
vTestData[2] = codeVersion;
vTestData[3] = cpuInfo;
vTestData[4] = benchDate;

vTestHeader[1] = "Solver Label";
vTestHeader[2] = "Solver Version";
vTestHeader[3] = "System Info";
vTestHeader[4] = "Test Date Time";

colIdx   = 5;
testIdx     = 1;

startTime = time();
for problemClass ∈ instances(ProblemClass)
    for iDim ∈ 1:size(mNumElements, 2)
        for iSim ∈ 1:numSimulations
            mP, vQ, mA, vL, vU = GenerateRandomQP(problemClass, mNumElements[Int(problemClass), iDim]; numConstraints = mNumConstraints[Int(problemClass), iDim]);
            # Some problems might change the actual data
            numElements     = size(mP, 1);
            numConstraints  = size(mA, 1);

            vXX = zeros(numElements);

            sBenchMark              = @benchmarkable $hSolverFun(vX, $mP, $vQ, $mA, $vL, $vU) setup = (vX = copy($vXX)); #<! We must interpolate the function as it is not in the global scope
            cBenchMark, convFlag    = BenchmarkTools.run_result(sBenchMark, samples = benchMarkSamples, evals = benchMarkEvals, seconds = benchMarkSeconds);

            vTestData[colIdx + 0] = min(cBenchMark.times...);
            vTestData[colIdx + 1] = cBenchMark.allocs;
            vTestData[colIdx + 2] = cBenchMark.memory;
            vTestData[colIdx + 3] = convFlag ≠ convNumItr;

            vTestHeader[colIdx + 0] = "Test $(@sprintf("%04d", testIdx)) Run Time";
            vTestHeader[colIdx + 1] = "Test $(@sprintf("%04d", testIdx)) # Allocations";
            vTestHeader[colIdx + 2] = "Test $(@sprintf("%04d", testIdx)) Allocations Size";
            vTestHeader[colIdx + 3] = "Test $(@sprintf("%04d", testIdx)) Convergence";

            println("Finished Test #$(@sprintf("%04d", testIdx)) Out of $(@sprintf("%04d", numTests)) Tests")

            global colIdx += 4;
            global testIdx += 1;
        end
    end
end
runTime = time() - startTime;

println("\nTotal Run Time: $(runTime) [Sec]");

mTestData = reshape(vTestData, 1, :);

if (isfile(benchmarkCsvFileName))
    csvFile = CSV.File(benchmarkCsvFileName);
    vCsvHeader = String.(propertynames(csvFile));
    if (vCsvHeader == vTestHeader)
        tTestData = Tables.table(mTestData); #<! No header
        CSV.write(benchmarkCsvFileName, tTestData; append = true);
    else
        error("The Header of the tests doesn't match the header of the CSV file");
    end
else
    tTestData = Tables.table(mTestData; header = vTestHeader);
    CSV.write(benchmarkCsvFileName, tTestData);
end

max_stars_repo_name: RoyiAvital/QuadraticProgramSolver

max_issues_repo_name: RoyiAvital/QuadraticProgramSolver

max_forks_repo_name: RoyiAvital/QuadraticProgramSolver


text:
# Top level Runner functions for gridsearch (optimization)
# for the Categorical task (Section 2.3.2)
# and for the robustness analysis for the Categorical task (Figure 9)

function runner_gridsearch_multinomial()
    run_gridsearch_multinomial_pf(20)
    run_gridsearch_multinomial_pf(20, nsteps = 10^6,
                                changeprobabilities = [0.001, 0.0001])
    run_gridsearch_multinomial_varsmile()
    run_gridsearch_multinomial_varsmile(nsteps = 10^6,
                                         changeprobabilities = [0.001, 0.0001])
end
function runner_gridrobustness_multinomial()
    run_gridsearch_multinomial_pf(20,
                                    nsteps = 10^5,
                                    stochasticities = [5., 0.14], # for robustness
                                    changeprobabilities = [10^i for i in range(-4, length=40, stop=-1)], # for robustness
                                    )
    run_gridsearch_multinomial_varsmile(
                                    nsteps = 10^5,
                                    stochasticities = [5., 0.14], # for robustness
                                    changeprobabilities = [10^i for i in range(-4, length=40, stop=-1)], # for robustness
                                    )
end

# -------------------------
function run_gridsearch_multinomial_leaky(;# nsteps = 10^6,
                                            nsteps = 10^5,
                                            # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                                            changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps)
                                            stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01], # Stays always same, except for robustness
                                            )
    learnertype = TLeakyIntegrator
    learnerparametervalues = collect(eps():0.01:0.9999999)
    foldername = "gridrobust_multi"
    #foldername = "grid_multi"
    filename = "leaky"
    run_gridsearch_multinomial(learnertype = learnertype,
                learnerparametervalues = learnerparametervalues,
                foldername = foldername,
                filename = filename,
                nsteps = nsteps,
                changeprobabilities = changeprobabilities,
                stochasticities = stochasticities,)
end
function run_gridsearch_multinomial_smileOriginal(;# nsteps = 10^6,
                                        nsteps = 10^5,
                                        # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                                        changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps)
                                        stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01], # Stays always same, except for robustness
                                        )
    learnertype = TSmileOriginal
    pc_log = [10^i for i in range(-4, length=40, stop=-1)]
    learnerparametervalues_linear = collect(eps():0.01:10.)
    learnerparametervalues = vcat(pc_log[1:end-1], learnerparametervalues_linear)
    foldername = "gridrobust_multi"
    #foldername = "grid_multi"
    filename = "smileOriginal_m"
    run_gridsearch_multinomial(learnertype = learnertype,
                learnerparametervalues = learnerparametervalues,
                foldername = foldername,
                filename = filename,
                nsteps = nsteps,
                changeprobabilities = changeprobabilities,
                stochasticities = stochasticities,)
end
function run_gridsearch_multinomial_varsmile(;# nsteps = 10^6,
                                            nsteps = 10^5,
                                            # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                                            changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps)
                                            stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01], # Stays always same, except for robustness
                                            )
    learnertype = TVarSmile
    pc_log = [10^i for i in range(-4, length=40, stop=-1)]
    learnerparametervalues_linear = collect(eps():0.01:10.)
    learnerparametervalues = vcat(pc_log[1:end-1], learnerparametervalues_linear)
    foldername = "gridrobust_multi"
    #foldername = "grid_multi"
    filename = "varsmile_m"
    run_gridsearch_multinomial(learnertype = learnertype,
                learnerparametervalues = learnerparametervalues,
                foldername = foldername,
                filename = filename,
                nsteps = nsteps,
                changeprobabilities = changeprobabilities,
                stochasticities = stochasticities,)
end
function run_gridsearch_multinomial_pf(Npar;# nsteps = 10^6,
                                        nsteps = 10^5,
                                        # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                                        changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps)
                                        stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01], # Stays always same, except for robustness
                                        )
    learnertype = TParticleFilter
    learnerparametervalues = [10^i for i in range(-4, length=40, stop=-1)]
    foldername = "gridrobust_multi"
    filename = "pf"*string(Npar)*"_pc"
    run_gridsearch_multinomial(learnertype = learnertype,
                learnerparametervalues = learnerparametervalues,
                foldername = foldername,
                filename = filename,
                Npar = Npar,
                nsteps = nsteps,
                changeprobabilities = changeprobabilities,
                stochasticities = stochasticities,)
end
function run_gridsearch_multinomial_tbayes(;# nsteps = 10^6,
                                        nsteps = 10^5,
                                        # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                                        changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps)
                                        stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01], # Stays always same, except for robustness
                                        )
    learnertype = TBayesFilter
    learnerparametervalues = [10^i for i in range(-4, length=40, stop=-1)]
    foldername = "gridrobust_multi"
    filename = "tbayes_pc"
    run_gridsearch_multinomial(learnertype = learnertype,
                learnerparametervalues = learnerparametervalues,
                foldername = foldername,
                filename = filename,
                nsteps = nsteps,
                changeprobabilities = changeprobabilities,
                stochasticities = stochasticities,)
end
function run_gridsearch_multinomial_tmpn(Npar;# nsteps = 10^6,
                                        nsteps = 10^5,
                                        # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                                        changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps)
                                        stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01], # Stays always same, except for robustness
                                        )
    learnertype = TMPN
    learnerparametervalues = [10^i for i in range(-4, length=40, stop=-1)]
    foldername = "gridrobust_multi"
    filename = "tsor"*string(Npar)*"_pc"
    run_gridsearch_multinomial(learnertype = learnertype,
                learnerparametervalues = learnerparametervalues,
                foldername = foldername,
                filename = filename,
                Npar = Npar,
                nsteps = nsteps,
                changeprobabilities = changeprobabilities,
                stochasticities = stochasticities,)
end
function run_gridsearch_multinomial_tsor(Npar;# nsteps = 10^6,
                                        nsteps = 10^5,
                                        # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                                        changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps)
                                        stochasticities = [5., 2., 1., 0.25, 0.14, 0.1, 0.01], # Stays always same, except for robustness
                                        )
    learnertype = TSOR
    learnerparametervalues = [10^i for i in range(-4, length=40, stop=-1)]
    foldername = "gridrobust_multi"
    filename = "tsorOriginal"*string(Npar)*"_pc"
    run_gridsearch_multinomial(learnertype = learnertype,
                learnerparametervalues = learnerparametervalues,
                foldername = foldername,
                filename = filename,
                Npar = Npar,
                nsteps = nsteps,
                changeprobabilities = changeprobabilities,
                stochasticities = stochasticities,)
end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the `EnhancedTestSet` type.

This set of unit tests checks the behavior of `EnhancedTestSet` with nested test sets.

Notes
-----
* For the unit tests in this files, failures and errors are expected.
"""

# --- Imports

# Standard library
using Test
using Test: DefaultTestSet

# External packages
using Suppressor

# Local modules
using TestTools.jltest

# --- Tests

# ------ Failing tests with diffs

# Nested inherited test set
output = strip(
    @capture_out begin
        @testset EnhancedTestSet "EnhancedTestSet: nested inherited EnhancedTestSet" begin
            @testset "Nested Inherited Test Set" begin
                @test [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]
            end
        end
    end
)

expected_prefix = strip("""
                  =====================================================
                  Nested Inherited Test Set: Test Failed
                    Expression: [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]

                    Diff:
                  [3, 5, 6, 1, (-)6, (+)9, 8]
                  """)
@test startswith(output, expected_prefix)

# Nested DefaultTestSet
output = strip(@capture_out begin
    @testset EnhancedTestSet "EnhancedTestSet: nested DefaultTestSet" begin
        @testset DefaultTestSet "DefaultTestSet Nested in EnhancedTestSet" begin
            @test [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]
        end
    end
end)

expected_prefix = strip(
    """
    DefaultTestSet Nested in EnhancedTestSet: Test Failed at $(@__FILE__):66
      Expression: [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]
       Evaluated: [3, 5, 6, 1, 6, 8] == [3, 5, 6, 1, 9, 8]
    Stacktrace:
    """
)
@test startswith(output, expected_prefix)

# --- Emit message about expected failures and errors

println()
@info "For $(basename(@__FILE__)), 2 failures and 0 error are expected."


text:
using Pkg
@static if isdefined(Base, :parsed_toml)
    parsed_toml(args...) = Base.parsed_toml(args...)
else
    parsed_toml(file) = Pkg.TOML.parsefile(file)
end

const project_names = ("JuliaProject.toml", "Project.toml")
const manifest_names = ("JuliaManifest.toml", "Manifest.toml")

# return nothing or the project file at env
function env_file(env::String, names = project_names)::Union{Nothing,String}
    if isdir(env)
        for proj in names
            project_file = joinpath(env, proj)
            safe_isfile(project_file) && return project_file
        end
        return nothing
    elseif basename(env) in names && safe_isfile(env)
        return env
    end
    return nothing
end

function is_project_folder_in_env(folder, env_manifest, server)
    project_file = Base.env_project_file(folder)
    project_file isa Bool && return false
    folder_proj = parsed_toml(project_file)
    manifest_pe = get(env_manifest, get(folder_proj, "uuid", ""), nothing)
    if manifest_pe === nothing
        return false
    elseif manifest_pe.path !== nothing
        return Base.Filesystem.samefile(manifest_pe.path, folder)
    elseif manifest_pe.tree_hash isa Base.SHA1
        if Base.Filesystem.samefile(abspath(server.depot_path, "packages", folder_proj.name, Base.version_slug(folder_proj.uuid, manifest_pe.tree_hash, 4)), folder)
            return true
        elseif Base.Filesystem.samefile(abspath(server.depot_path, "packages", folder_proj.name, Base.version_slug(folder_proj.uuid, manifest_pe.tree_hash)), folder)
            return true
        end
    end
    return false
end

function get_env_for_root(doc::Document, server::LanguageServerInstance)
    env_proj_file = env_file(server.env_path, project_names)
    env_manifest_file = env_file(server.env_path, manifest_names)

    (safe_isfile(env_proj_file) && safe_isfile(env_manifest_file)) || return

    # Find which workspace folder the doc is in.
    parent_workspaceFolders = sort(filter(f->startswith(doc._path, f), collect(server.workspaceFolders)), by = length, rev = true)

    isempty(parent_workspaceFolders) && return
    # arbitrarily pick one
    parent_workspaceFolder = first(parent_workspaceFolders)

    project_env = env_file(parent_workspaceFolder, project_names)
    try
        if safe_isfile(project_env)
            folder_proj = parsed_toml(project_env)

            # We point to all caches as, though a package may not be directly available (e.g as
            # a dependency) it may still be accessible as imported by one of the direct dependencies.
            symbols = server.global_env.symbols

            # Will want to limit this to only get extended methods from the dependency tree of
            # the project (e.g. using `complete_dep_tree` below)
            extended_methods = server.global_env.extended_methods

            # This is the list of packages that are directly available
            project_deps = Symbol.(collect(keys(get(folder_proj, "deps", []))))
            if isdir(joinpath(parent_workspaceFolder, "test")) && startswith(doc._path, joinpath(parent_workspaceFolder, "test"))
                # We're in the test folder, add the project iteself to the deps
                # This should actually point to the live code (e.g. the relevant EXPR or Scope)?
                haskey(folder_proj, "name") && push!(project_deps, Symbol(folder_proj["name"]))
                for extra in keys(get(folder_proj, "extras", []))
                    if Symbol(extra) in keys(symbols)
                        push!(project_deps, Symbol(extra))
                    end
                end
            end

            StaticLint.ExternalEnv(symbols, extended_methods, project_deps)
        end
    catch err
        # The specified env is faulty (incorrect format, missing entries, ...)
        # We can't do anything about that except treating it the same as a non-existing env,
        # but showing a warning might be useful.
        msg = "The Julia environment at `$project_env` is invalid. Using the global environment instead."
        if server.jr_endpoint !== nothing
            JSONRPC.send(server.jr_endpoint, window_showMessage_notification_type, ShowMessageParams(
                MessageTypes.Warning,
                msg
            ))
        end
        @error msg exception=(err, catch_backtrace())
    end
end

function complete_dep_tree(uuid, env_manifest, alldeps = Dict{Base.UUID,Pkg.Types.PackageEntry}())
    haskey(alldeps, uuid) && return alldeps
    alldeps[uuid] = env_manifest[uuid]
    for dep_uuid in values(alldeps[uuid].deps)
        complete_dep_tree(dep_uuid, env_manifest, alldeps)
    end
    alldeps
end


text:
export GroundStateProblem

"""
    GroundStateProblem([T=STD_REAL_PREC], ham, operators=true, variance=false)

Returns the problem targeting of the ground state energy. If `variance=false`
this is done by minimizing the energy, otherwise it's done by minimizing the
variance.

If `operators=true` a memory efficient representation of the hamiltonian is used,
resulting in less memory consuption but higher CPU usage. This is needed for lattices
bigger than a certain threshold.
To use this feature, the hamiltonian must be provided as a `GraphOperator` object,
from `QuantumLattices` package.
"""
GroundStateProblem(args...; kwargs...) = GroundStateProblem(STD_REAL_PREC, args...; kwargs...)
GroundStateProblem(T::Type{<:Number}, args...; kwargs...) = _build_groundstate_problem(T, args...; kwargs...)

# If Called with QuantumOptics' DataOperators, we can't convert them to
# KLocalOperators, so set `operators=false` and warn the user
function GroundStateProblem(T::Type{<:Number}, hamiltonian::DataOperator; operators=false, kwargs...)
    if operators
        @warn """Converting QuantumOptics' DataOperators to memory-efficient Local Operators is not possible.
        Setting `operators=false`.
        To suppress this warning, don't request `operators=true`.
        """
        operators = false
    end
    return _build_groundstate_problem(T, hamiltonian, operators=operators; kwargs...)
end

# The actual implementation
function _build_groundstate_problem(T::Type{<:Number}, hamiltonian; operators=true, variance=false)
    base = basis(hamiltonian)
    T    = Complex{real(T)}

    if operators
        ham = to_linear_operator(hamiltonian, T)
    else
        ham = convert(AbstractMatrix{T},
                      data(SparseOperator(hamiltonian)))
    end

    if !variance
        return HamiltonianGSEnergyProblem(basis(hamiltonian), ham, 0.0)
    else
        throw("Variance minimization problem not implemented yet")
    end
end


text:
# Adapt method for DynamicGrids objects
function Adapt.adapt_structure(to, x::AbstractSimData)
    @set! x.grids = map(g -> Adapt.adapt(to, g), x.grids)
    @set! x.extent = Adapt.adapt(to, x.extent)
    return x
end

function Adapt.adapt_structure(to, x::GridData)
    @set! x.source = Adapt.adapt(to, x.source)
    @set! x.source = Adapt.adapt(to, x.source)
    @set! x.mask = Adapt.adapt(to, x.mask)
    @set! x.dest = Adapt.adapt(to, x.dest)
    @set! x.optdata = Adapt.adapt(to, x.optdata)
    return x
end

function Adapt.adapt_structure(to, x::AbstractExtent)
    @set! x.init = _adapt_x(to, init(x))
    @set! x.mask = _adapt_x(to, mask(x))
    @set! x.aux = _adapt_x(to, aux(x))
    return x
end

_adapt_x(to, A::AbstractArray) = Adapt.adapt(to, A)
_adapt_x(to, nt::NamedTuple) = map(A -> Adapt.adapt(to, A), nt)
_adapt_x(to, nt::Nothing) = nothing

# Adapt output frames to GPU
# TODO: this may be incorrect use of Adapt.jl, as the Output
# object is not entirely adopted for GPU use, the CuArray
# frames are still held in a regular Array.
function Adapt.adapt_structure(to, o::Output)
    frames = map(o.frames) do f
        _adapt_x(to, f)
    end
    @set! o.extent = adapt(to, o.extent)
    @set! o.frames = frames
    return o
end


text:
include("solver_preamble.jl")
import SDPA
factory = with_optimizer(SDPA.Optimizer)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5, query=false)
@testset "Linear" begin
    # With `dsos_horn_test`, the termination status is
    # `INFEASIBLE_OR_UNBOUNDED` instead of `INFEASIBLE`.
    # With `dsos_concave_then_convex_cubic`, we have
    # cholesky miss condition :: not positive definite :: info = 32 :: line 785 in sdpa_linear.cpp
    # There are some possibilities. :: line 786 in sdpa_linear.cpp
    # 1. SDPA terminates due to inaccuracy of numerical error :: line 787 in sdpa_linear.cpp
    # 2. The input problem may not have (any) interior-points :: line 788 in sdpa_linear.cpp
    # 3. Input matrices are linearly dependent :: line 789 in sdpa_linear.cpp
    Tests.linear_test(factory, config,
                      ["dsos_horn", "dsos_concave_then_convex_cubic"])
end
@testset "SDP" begin
    # With `sos_concave_then_convex_cubic`, we have the same error than with
    # `dsos_concave_then_convex_cubic`.
    Tests.sd_test(factory, config,
                  ["sos_concave_then_convex_cubic",
                   "quadratic_infeasible_lyapunov_switched_system",
                   "quartic_infeasible_lyapunov_switched_system"])
end


text:
incx < 0 && (ix = ix+(n-1)*abs(incx))
incy < 0 && (iy = iy+(n-1)*abs(incy))
incz < 0 && (iz = iz+(n-1)*abs(incz))
@inbounds for i = 0:n-1
    x[ix+i*incx] = y[iy+i*incy] + z[iz+i*incz]*a
end


text:
"""
    test_scalar(f, z; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), fkwargs=NamedTuple(), check_inferred=true, kwargs...)

Given a function `f` with scalar input and scalar output, perform finite differencing checks,
at input point `z` to confirm that there are correct `frule` and `rrule`s provided.

# Arguments
- `f`: function for which the `frule` and `rrule` should be tested.
- `z`: input at which to evaluate `f` (should generally be set to an arbitrary point in the domain).

# Keyword Arguments
- `fdm`: the finite differencing method to use.
- `fkwargs` are passed to `f` as keyword arguments.
- If `check_inferred=true`, then the inferrability (type-stability) of the `frule` and `rrule` are checked.
- All remaining keyword arguments are passed to `isapprox`.
"""
function test_scalar(f, z; rtol=1e-9, atol=1e-9, fdm=_fdm, fkwargs=NamedTuple(), check_inferred=true, kwargs...)
    # To simplify some of the calls we make later lets group the kwargs for reuse
    rule_test_kwargs = (; rtol=rtol, atol=atol, fdm=fdm, fkwargs=fkwargs, check_inferred=check_inferred, kwargs...)
    isapprox_kwargs = (; rtol=rtol, atol=atol, kwargs...)

    @testset "test_scalar: $f at $z" begin
        # z = x + im * y
        # Ω = u(x, y) + im * v(x, y)
        Ω = f(z; fkwargs...)

        # test jacobian using forward mode
        Δx = one(z)
        @testset "with tangent $Δx" begin
            # check ∂u_∂x and (if Ω is complex) ∂v_∂x via forward mode
            test_frule(f, z ⊢ Δx; rule_test_kwargs...)
            if z isa Complex
                # check that same tangent is produced for tangent 1.0 and 1.0 + 0.0im
                ḟ = rand_tangent(f)
                _, real_tangent = frule((ḟ, real(Δx)), f, z; fkwargs...)
                _, embedded_tangent = frule((ḟ, Δx), f, z; fkwargs...)
                test_approx(real_tangent, embedded_tangent; isapprox_kwargs...)
            end
        end
        if z isa Complex
            Δy = one(z) * im
            @testset "with tangent $Δy" begin
                # check ∂u_∂y and (if Ω is complex) ∂v_∂y via forward mode
                test_frule(f, z ⊢ Δy; rule_test_kwargs...)
            end
        end

        # test jacobian transpose using reverse mode
        Δu = one(Ω)
        @testset "with cotangent $Δu" begin
            # check ∂u_∂x and (if z is complex) ∂u_∂y via reverse mode
            test_rrule(f, z ⊢ Δx; output_tangent=Δu, rule_test_kwargs...)
            if Ω isa Complex
                # check that same cotangent is produced for cotangent 1.0 and 1.0 + 0.0im
                _, back = rrule(f, z)
                _, real_cotangent = back(real(Δu))
                _, embedded_cotangent = back(Δu)
                test_approx(real_cotangent, embedded_cotangent; isapprox_kwargs...)
            end
        end
        if Ω isa Complex
            Δv = one(Ω) * im
            @testset "with cotangent $Δv" begin
                # check ∂v_∂x and (if z is complex) ∂v_∂y via reverse mode
                test_rrule(f, z ⊢ Δx; output_tangent=Δv, rule_test_kwargs...)
            end
        end
    end  # top-level testset
end

"""
    test_frule([config::RuleConfig,] f, args..; kwargs...)

# Arguments
- `config`: defaults to `ChainRulesTestUtils.ADviaRuleConfig`.
- `f`: function for which the `frule` should be tested. Its tangent can be provided using `f ⊢ ḟ`.
  (You can enter `⊢` via `\\vdash` + tab in the Julia REPL and supporting editors.)
- `args...`: either the primal args `x`, or primals and their tangents: `x ⊢ ẋ`
   - `x`: input at which to evaluate `f` (should generally be set to an arbitrary point in the domain).
   - `ẋ`: differential w.r.t. `x`; will be generated automatically if not provided.
   Non-differentiable arguments, such as indices, should have `ẋ` set as `NoTangent()`.

# Keyword Arguments
   - `output_tangent`: tangent against which to test accumulation of derivatives.
     Should be a differential for the output of `f`. Is set automatically if not provided.
   - `fdm::FiniteDifferenceMethod`: the finite differencing method to use.
   - `frule_f=frule`: function with an `frule`-like API that is tested (defaults to
     `frule`). Used for testing gradients from AD systems.
   - If `check_inferred=true`, then the inferrability (type-stability) of the `frule` is checked,
     as long as `f` is itself inferrable.
   - `fkwargs` are passed to `f` as keyword arguments.
   - All remaining keyword arguments are passed to `isapprox`.
"""
function test_frule(args...; kwargs...)
    config = ChainRulesTestUtils.ADviaRuleConfig()
    test_frule(config, args...; kwargs...)
end

function test_frule(
    config::RuleConfig,
    f,
    args...;
    output_tangent=Auto(),
    fdm=_fdm,
    frule_f=ChainRulesCore.frule,
    check_inferred::Bool=true,
    fkwargs::NamedTuple=NamedTuple(),
    rtol::Real=1e-9,
    atol::Real=1e-9,
    kwargs...,
)
    # To simplify some of the calls we make later lets group the kwargs for reuse
    isapprox_kwargs = (; rtol=rtol, atol=atol, kwargs...)

    # and define a helper closure
    call_on_copy(f, xs...) = deepcopy(f)(deepcopy(xs)...; deepcopy(fkwargs)...)

    @testset "test_frule: $f on $(_string_typeof(args))" begin

        primals_and_tangents = auto_primal_and_tangent.((f, args...))
        primals = primal.(primals_and_tangents)
        tangents = tangent.(primals_and_tangents)

        if check_inferred && _is_inferrable(deepcopy(primals)...; deepcopy(fkwargs)...)
            _test_inferred(frule_f, deepcopy(config), deepcopy(tangents), deepcopy(primals)...; deepcopy(fkwargs)...)
        end

        res = call_on_copy(frule_f, config, tangents, primals...)
        res === nothing && throw(MethodError(frule_f, typeof(primals)))
        @test_msg "The frule should return (y, ∂y), not $res." res isa Tuple{Any,Any}
        Ω_ad, dΩ_ad = res
        Ω = call_on_copy(primals...)
        test_approx(Ω_ad, Ω; isapprox_kwargs...)

        # Correctness testing via finite differencing.
        is_ignored = isa.(tangents, NoTangent)
        dΩ_fd = _make_jvp_call(fdm, call_on_copy, Ω, primals, tangents, is_ignored)
        test_approx(dΩ_ad, dΩ_fd; isapprox_kwargs...)

        acc = output_tangent isa Auto ? rand_tangent(Ω) : output_tangent
        _test_add!!_behaviour(acc, dΩ_ad; isapprox_kwargs...)
    end  # top-level testset
end

"""
    test_rrule([config::RuleConfig,] f, args...; kwargs...)

# Arguments
- `config`: defaults to `ChainRulesTestUtils.ADviaRuleConfig`.
- `f`: function for which the `rrule` should be tested. Its tangent can be provided using `f ⊢ f̄`.
  (You can enter `⊢` via `\\vdash` + tab in the Julia REPL and supporting editors.)
- `args...`: either the primal args `x`, or primals and their tangents: `x ⊢ x̄`
    - `x`: input at which to evaluate `f` (should generally be set to an arbitrary point in the domain).
    - `x̄`: currently accumulated cotangent; will be generated automatically if not provided.
    Non-differentiable arguments, such as indices, should have `x̄` set as `NoTangent()`.

# Keyword Arguments
 - `output_tangent`: the seed to propagate backward for testing (technically a cotangent).
   should be a differential for the output of `f`. Is set automatically if not provided.
 - `check_thunked_output_tangent=true`: also checks that passing a thunked version of the
    output tangent to the pullback returns the same result.
 - `fdm::FiniteDifferenceMethod`: the finite differencing method to use.
 - `rrule_f=rrule`: function with an `rrule`-like API that is tested (defaults to `rrule`).
   Used for testing gradients from AD systems.
 - If `check_inferred=true`, then the inferrability (type-stability) of the `rrule` is checked
   — if `f` is itself inferrable — along with the inferrability of the pullback it returns.
 - `fkwargs` are passed to `f` as keyword arguments.
 - All remaining keyword arguments are passed to `isapprox`.
"""
function test_rrule(args...; kwargs...)
    config = ChainRulesTestUtils.ADviaRuleConfig()
    test_rrule(config, args...; kwargs...)
end

function test_rrule(
    config::RuleConfig,
    f,
    args...;
    output_tangent=Auto(),
    check_thunked_output_tangent=true,
    fdm=_fdm,
    rrule_f=ChainRulesCore.rrule,
    check_inferred::Bool=true,
    fkwargs::NamedTuple=NamedTuple(),
    rtol::Real=1e-9,
    atol::Real=1e-9,
    kwargs...,
)
    # To simplify some of the calls we make later lets group the kwargs for reuse
    isapprox_kwargs = (; rtol=rtol, atol=atol, kwargs...)

    # and define helper closure over fkwargs
    call(f, xs...) = f(xs...; fkwargs...)

    @testset "test_rrule: $f on $(_string_typeof(args))" begin

        # Check correctness of evaluation.
        primals_and_tangents = auto_primal_and_tangent.((f, args...))
        primals = primal.(primals_and_tangents)
        accum_cotangents = tangent.(primals_and_tangents)

        if check_inferred && _is_inferrable(primals...; fkwargs...)
            _test_inferred(rrule_f, config, primals...; fkwargs...)
        end
        res = rrule_f(config, primals...; fkwargs...)
        res === nothing && throw(MethodError(rrule_f, typeof(primals)))
        y_ad, pullback = res
        y = call(primals...)
        test_approx(y_ad, y; isapprox_kwargs...)  # make sure primal is correct

        ȳ = output_tangent isa Auto ? rand_tangent(y) : output_tangent

        check_inferred && _test_inferred(pullback, ȳ)
        ad_cotangents = pullback(ȳ)
        @test_msg(
            "The pullback must return a Tuple (∂self, ∂args...)",
            ad_cotangents isa Tuple
        )
        @test_msg(
            "The pullback should return 1 cotangent for the primal and each primal input.",
            length(ad_cotangents) == length(primals)
        )

        # Correctness testing via finite differencing.
        is_ignored = isa.(accum_cotangents, NoTangent)
        fd_cotangents = _make_j′vp_call(fdm, call, ȳ, primals, is_ignored)
        foreach(accum_cotangents, ad_cotangents, fd_cotangents) do args...
            _test_cotangent(args...; check_inferred=check_inferred, isapprox_kwargs...)
        end

        if check_thunked_output_tangent
            test_approx(ad_cotangents, pullback(@thunk(ȳ)), "pulling back a thunk:")
            check_inferred && _test_inferred(pullback, @thunk(ȳ))
        end
    end  # top-level testset
end

"""
    @maybe_inferred [Type] f(...)

Like `@inferred`, but does not check the return type if tests are run as part of PkgEval or
if the environment variable `CHAINRULES_TEST_INFERRED` is set to `false`.
"""
macro maybe_inferred(ex...)
    inferred = Expr(:macrocall, GlobalRef(Test, Symbol("@inferred")), __source__, ex...)
    return :(TEST_INFERRED[] ? $(esc(inferred)) : $(esc(last(ex))))
end

"""
    _test_inferred(f, args...; kwargs...)

Simple wrapper for [`@maybe_inferred f(args...: kwargs...)`](@ref `@maybe_inferred`), avoiding the type-instability in not
knowing how many `kwargs` there are.
"""
function _test_inferred(f, args...; kwargs...)
    if isempty(kwargs)
        @maybe_inferred f(args...)
    else
        @maybe_inferred f(args...; kwargs...)
    end
end

"""
    _is_inferrable(f, args...; kwargs...) -> Bool

Return whether the return type of `f(args...; kwargs...)` is inferrable.
"""
function _is_inferrable(f, args...; kwargs...)
    try
        _test_inferred(f, args...; kwargs...)
        return true
    catch ErrorException
        return false
    end
end

"""
    _test_cotangent(accum_cotangent, ad_cotangent, fd_cotangent; kwargs...)

Check if the cotangent `ad_cotangent` from `rrule` is consistent with `accum_tangent` and
approximately equal to the cotangent `fd_cotangent` obtained with finite differencing.

If `accum_cotangent` is `NoTangent()`, i.e., the argument was marked as non-differentiable,
`ad_cotangent` and `fd_cotangent` should be `NoTangent()` as well.

# Keyword arguments
- If `check_inferred=true` (the default) and `ad_cotangent` is a thunk, then it is checked if
  its content can be inferred.
- All remaining keyword arguments are passed to `isapprox`.
"""
function _test_cotangent(
    accum_cotangent,
    ad_cotangent,
    fd_cotangent;
    check_inferred=true,
    kwargs...,
)
    ad_cotangent isa AbstractThunk && check_inferred && _test_inferred(unthunk, ad_cotangent)

    # The main test of the actual derivative being correct:
    test_approx(ad_cotangent, fd_cotangent; kwargs...)
    _test_add!!_behaviour(accum_cotangent, ad_cotangent; kwargs...)
end

# we marked the argument as non-differentiable
function _test_cotangent(::NoTangent, ad_cotangent, ::NoTangent; kwargs...)
    @test ad_cotangent isa NoTangent
end
function _test_cotangent(::NoTangent, ::ZeroTangent, ::NoTangent; kwargs...)
    error(
        "The pullback in the rrule should use NoTangent()" *
        " rather than ZeroTangent() for non-perturbable arguments."
    )
end
function _test_cotangent(
    ::NoTangent,
    ad_cotangent::ChainRulesCore.NotImplemented,
    ::NoTangent;
    kwargs...,
)
    # this situation can occur if a cotangent is not implemented and
    # the default `rand_tangent` is `NoTangent`: e.g. due to having no fields
    # the `@test_broken` below should tell them that there is an easy implementation for
    # this case of `NoTangent()` (`@test_broken false` would be less useful!)
    # https://github.com/JuliaDiff/ChainRulesTestUtils.jl/issues/217
    @test_broken ad_cotangent isa NoTangent
end
function _test_cotangent(::NoTangent, ad_cotangent, fd_cotangent; kwargs...)
    error("cotangent obtained with finite differencing has to be NoTangent()")
end


text:
struct SimpleLogger{T}
    data::T
    f_map::Dict{Symbol, Function}
    function SimpleLogger(names, dtypes, f_map)
        nt = NamedTuple{names}([Vector{dt}() for dt ∈ dtypes])
        new{typeof(nt)}(nt,
                        f_map)
    end
end

get_data(sl::SimpleLogger) = sl.data
Base.getindex(sl::SimpleLogger, idx::Symbol) = sl.data[idx]

function (sl::SimpleLogger)(args...; kwargs...)
    for k ∈ keys(sl.f_map)
        push!(getindex(sl.data, k), sl.f_map[k](args...; kwargs...))
    end
end

mutable struct UpdateStateAnalysis{T}
    data::T
    f_map::Dict{Symbol, Function}
    function UpdateStateAnalysis(initial, f_map)
        nt = initial
        new{typeof(nt)}(
            nt,
            f_map
        )
    end
end

function (usa::UpdateStateAnalysis)(us::UpdateState)
    usa.data = typeof(usa.data)((usa.f_map[n](getindex(usa.data, n), us)) for n ∈ keys(usa.data))
end

Base.getindex(usa::UpdateStateAnalysis, idx::Symbol) = getindex(usa.data, idx)

function l1_grad(s, us::UpdateState)
    ns = zero(s)
    for k ∈ us.grads.params
        if !(us.grads[k] isa Nothing)
            ns += sum(abs.(us.grads[k]))
        end
    end
    ns
end

function l2_grad(s, us::UpdateState)
    ns = zero(s)
    for (k, v) ∈ us.grads
        ns += dot(v, v)
    end
    ns
end


text:
#**************************************************************************************
# test_Input_Interface.jl
# ======================= part of the GeoEfficiency.jl package.
#
# test all the input to the package from either the console or the csv files.
#
#**************************************************************************************

@testset "Input_Interface" begin

print("\n\t"); @info("test `setSrcToPoint` & `typeofSrc`...")
	@testset "setSrcToPoint" begin
		@test G.srcType === G.srcUnknown  	# the initial value
		@test typeofSrc() === G.srcUnknown  	# the initial value
		@test setSrcToPoint() === false      # not defined, set to not point

		@test setSrcToPoint(false) === false
		@test G.srcType == G.srcNotPoint
		@test typeofSrc() === G.srcNotPoint
		@test setSrcToPoint() === false
		@test setSrcToPoint("\n Is it a point source {Y|n} ? ") === false

		@test setSrcToPoint(true) === true
		@test G.srcType === G.srcPoint
		@test typeofSrc() === G.srcPoint
		@test setSrcToPoint() === true
		@test setSrcToPoint("\n Is it a point source {Y|n} ? ") === true

		@test typeofSrc(-5) === G.srcUnknown
		@test typeofSrc(-1) === G.srcUnknown
		@test setSrcToPoint() === false
		#@test setSrcToPoint("\n Is it a point source {Y|n} ? ") === false # require input

		@test typeofSrc(0) === G.srcPoint
		@test setSrcToPoint() === true
		@test setSrcToPoint("\n Is it a point source {Y|n} ? ") === true
		@test setSrcToPoint(false) === false
		@test setSrcToPoint(true)  === true

		@test typeofSrc(1) === G.srcLine
		@test setSrcToPoint() === false
		@test setSrcToPoint("\n Is it a point source {Y|n} ? ") === false
		@test setSrcToPoint(false) === false
		@test setSrcToPoint(true)  === true

		@test typeofSrc(2) === G.srcDisk
		@test setSrcToPoint() === false
		@test setSrcToPoint("\n Is it a point source {Y|n} ? ") === false
		@test setSrcToPoint(false) === false
		@test setSrcToPoint(true)  === true

		@test typeofSrc(3) === G.srcVolume
		@test setSrcToPoint() === false
		@test setSrcToPoint("\n Is it a point source {Y|n} ? ") === false
		@test setSrcToPoint(false) === false
		@test setSrcToPoint(true)  === true

		@test typeofSrc(4) === G.srcNotPoint
		@test setSrcToPoint() === false
		@test setSrcToPoint("\n Is it a point source {Y|n} ? ") === false
		@test setSrcToPoint(false) === false
		@test setSrcToPoint(true)  === true

		@test typeofSrc(5) === G.srcNotPoint
		@test setSrcToPoint() === false
		@test setSrcToPoint(false) === false
		@test setSrcToPoint(true)  === true
	end #testset

print("\n\t"); @info("test `input`...")
	@testset "input" begin
		if true #!isapple()
			write( stdin.buffer," anbfyiQERFC \n")
			@test G.input() == " anbfyiQERFC "
		else
			@test_throws 	ErrorException	write(stdin.buffer," anbfyiQERFC \n")
		end #if
	end #testset_input

print("\n\t"); @info("test `getfloat`...")
	@testset "getfloat" begin
print("\n\t"); @info("test `getfloat` with different ways to input numbers...")
		@test   0.0     ==  getfloat("\njust press return: ",value="0.0")
		@test   1.0     ==  getfloat("\ninput 1, then press return: ",value="1.0")
		@test   1.0     ==  getfloat("\ninput 1.0, then press return: ",value="1.0")
		@test   2000.0  ==  getfloat("\ninput '2e3', then press return: ",value="2e3")
		@test   0.034   ==  getfloat("\ninput '3.4e-2', then press return: ",value="3.4e-2")
		@test   isa( getfloat("\ntry to input any string, only valid number should accepted: ",value="1*0im"), Float64)

print("\n\t"); @info("test `getfloat` with mathematical expressions...")
		@test   0.5           		==  getfloat("\ninput 1/2, then press return: ",value="1/2")
		@test   0.75          		==  getfloat("\ninput 3//4, then press return: ",value="3//4")
		@test   MathConstants.pi/2 	≈   getfloat("\ninput 'pi/2', then press return: ",value="pi/2")
		@test   MathConstants.e     ≈   getfloat("\ninput 'e', then press return: ",value="e")
		@test   MathConstants.e^3   ≈   getfloat("\ninput 'e^3', then press return: ",value="e^3")
		@test   Base.sin(0.1) 		≈   getfloat("\ninput 'sin(0.1)', then press return: ",value="sin(0.1)")
		if true #!isapple()
			for i = 0:5
				write(stdin.buffer,"1.2+2im\n"^i * "3\n")
				@test   3.0 == getfloat("\nthe first time input '1.2+2im': ")
			end # for
			write(stdin.buffer,"5\n" * "3\n")
			@test   3.0 == getfloat("\ninput 1/2, then press return: ", 0.0, 4.0)
			write(stdin.buffer,"-1\n" * "3\n")
			@test   3.0 == getfloat("\ninput 1/2, then press return: ", 0.0, 4.0)
			write(stdin.buffer,"1.2f\n" * "3\n")
			@test   3.0 == getfloat("\nthe first time input '1.2f': ", 0.0, 4.0)
		else
			@test_throws 	ErrorException	write(stdin.buffer,"5\n" * "3\n")
		end #if
  end # testset

print("\n\t"); @info("test `reading from CSV`...")
    @testset "reading from CSV" begin
		detector_info_array = [5 0 0 0; 5 10 0 0; 5 10 2 0; 5 10 2 5]
        detectors = [Detector(5, 0, 0, 0), Detector(5, 10, 0, 0), Detector(5, 10, 2, 0), Detector(5, 10, 2, 5)]

		datadirectory = joinpath(homedir(), "GeoEfficiency", "temptemp"); isdir(datadirectory) || mkdir(datadirectory)

		detectorfile = joinpath(datadirectory, "_Detector_test.csv")
		hightfile    = joinpath(datadirectory, "_hight_test.csv")
		Rhosfile     = "Rhosfile "
		Radiifile    = "Radiifile"
		Lengthsfile  = "Lengthsfile"
		setSrcToPoint(true) == true

	print("\n\t\t"); @info("Detectors write and read  - input type{Int}")
		@test  G.writecsv_head(detectorfile, detector_info_array, ["CryRadius"	 "CryLength" "HoleRadius" "HoleDepth"])  ==  nothing
		@test  G.detector_info_from_csvFile("_Detector_test.csv", datadirectory) == sort(detectors)

	print("\n\t\t"); @info("write and read  - input type{Int}")
		@test  G.writecsv_head(hightfile, [0, 1, 2, 3, 4, 5, 10, 15, 20,], ["SrcHight"])  ==  nothing
		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) == [0, 1, 2, 3, 4, 5, 10, 15, 20,]

	print("\n\t\t"); @info("READ_BATCH_INFO")
		batch_info = G.read_batch_info(datadirectory, detectorfile, hightfile, Rhosfile, Radiifile, Lengthsfile)
		@test  batch_info[1] == sort(detectors)
		@test  batch_info[2] == [0, 1, 2, 3, 4, 5, 10, 15, 20,]
		@test  batch_info[3] == [0.0]
		@test  batch_info[4] == [0.0]
		@test  batch_info[5] == [0.0]
		@test  batch_info[6] == ( G.srcType === G.srcPoint)
		@test  G.read_batch_info(datadirectory, detectorfile, hightfile, Rhosfile, Radiifile, Lengthsfile) == (detectors |> sort, [0.0, 1, 2, 3, 4, 5, 10, 15, 20,], [0.0], [0.0], [0.0], G.srcType === G.srcPoint)

	print("\n\t\t"); @info("rewrite, read and sort  - input type{Int}")
		@test  G.writecsv_head(hightfile, [3, 20, 4, 0, 1, 2, 5, 10, 15,], ["SrcHight"])  ==  nothing
		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) == [0, 1, 2, 3, 4, 5, 10, 15, 20,]

	print("\n\t\t"); @info("rewrite, read and sort - input type{Rational-treated as any}")
		@test  G.writecsv_head(hightfile, [3//2, 20, 4, 0, 1, 2, 5, 10, 15,], ["SrcHight"])  ==  nothing
		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) == [0.0]

	print("\n\t\t"); @info("rewrite, read and sort - input type{Float64}")
		@test  G.writecsv_head(hightfile, [3.0, 20, 4, 0, 1, 2, 5, 10, 15,], ["SrcHight"])  ==  nothing
		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) == [0, 1, 2, 3, 4, 5, 10, 15, 20,]

	#print("\n\t\t"); @info("rewrite, read and sort - input type{Irrational....}")
	#	@test  G.writecsv_head(hightfile, [pi, e, pi + e, 1, 2, 5, 10, 15,], ["SrcHight"])  ==  nothing
	#	@test G.read_from_csvFile("_hight_test.csv", datadirectory) == [0]

	print("\n\t\t"); @info("invalid data type {Unionall}")
		@test  G.writecsv_head(hightfile, ["3.0", 20, 4, 0, 1, 2, 5, 10, 15,], ["SrcHight"])  ==  nothing
		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) ==  [0, 1, 2, 3, 4, 5, 10, 15, 20,]

	print("\n\t\t"); @info("invalid data type {String}")
		@test  G.writecsv_head(hightfile, ["pi", "20", "4", "0", "1", "2", "5", "10", "15",], ["SrcHight"])  ==  nothing
		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) == [0.0]

	print("\n\t\t"); @info("invalid data type {Complex}")
		@test  G.writecsv_head(hightfile, [3.0+0.0im, 20, 4, 0, 1, 2, 5, 10, 15,], ["SrcHight"])  ==  nothing
		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) == [0.0]

	print("\n\t\t"); @info("missing file - `Hights.csv`")
		rm(hightfile, recursive=true)

		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) == [0.0]
		if true #!isapple()
			setSrcToPoint(false)
			write(stdin.buffer,
			"1\n" * "0\n" * #=axial point=#
			"2\n" * "3\n" #=SrcRadius SrcHeight=#)
			@test  G.read_batch_info(datadirectory, detectorfile, hightfile, Rhosfile, Radiifile, Lengthsfile) == (detectors |> sort, [1.0],	[0.0], [2.0], [3.0], false)
		end #if

	print("\n\t\t"); @info("missing file - `Rhos.csv`")
		@test  G.writecsv_head(hightfile, [4.0,], ["SrcHight"])  ==  nothing
		@test  G.read_from_csvFile("_hight_test.csv", datadirectory) == [4.0]
		if true #!isapple()
			setSrcToPoint(false)
			write(stdin.buffer,
			"1\n" * "0\n" * #=axial point=#
			"2\n" * "3\n" #=SrcRadius SrcHeight=#)
			@test  G.read_batch_info(datadirectory, detectorfile, hightfile, Rhosfile, Radiifile, Lengthsfile) == (detectors |> sort, [1.0],	[0.0], [2.0], [3.0], false)
		end #if

	print("\n\t\t"); @info("Detectors - missing file\n")
		rm(datadirectory, recursive=true)
		@test_throws Union{ArgumentError, SystemError}  G.detector_info_from_csvFile("_Detector_test.csv", datadirectory)   # the Union{ArgumentError, SystemError} is used for compatibility in both 0.6 and 0.7-dev

	rm(datadirectory, force=true, recursive=true)

	@test_throws Union{ArgumentError, SystemError}  G.detector_info_from_csvFile("jskfdsiu.uty","fghpweuh.uty")

	try
		@test eltype(G.detector_info_from_csvFile()) == G.Detector

		if  [0.0] != G.read_from_csvFile(G.srcHeights, G.datadir)
			setSrcToPoint(true);
			@test G.read_batch_info()[end]=== true

			if [0.0] != G.read_from_csvFile(G.srcRadii, G.datadir)
				setSrcToPoint(false);
				@test G.read_batch_info()[end]=== false
			end
		end
	end

	end # testset

print("\n\t"); @info("test `getDetectors`...")
    @testset "getDetectors" begin
		detector_info_array = [5 0 0 0; 5 10 0 0; 5 10 2 0; 5 10 2 5]
        detectors = [Detector(5, 0, 0, 0), Detector(5, 10, 0, 0), Detector(5, 10, 2, 0), Detector(5, 10, 2, 5)]
	    detectors = detectors |> sort

		@test getDetectors(detector_info_array) == detectors
		for det = detectors
			@test typeof(det) != Detector
		end
		@test eltype(detectors) != CylDetector
		@test eltype(detectors) == Detector
		det1, det2, det3, det4 = detectors
		@test det1 <= det2 <= det3 <= det4

		detector_info_array = [5 0; 10 0; 15 0; 20 0]
		detectors = getDetectors(detector_info_array)
		for det = detectors
			@test typeof(det) == CylDetector
			@test typeof(det) != Detector
		end
		@test eltype(detectors) != CylDetector
		@test eltype(detectors) == Detector
		det1, det2, det3, det4 = detectors

		detector_info_array = [5 1; 10 1; 15 1; 20 1]
		detectors = getDetectors(detector_info_array)
		for det = detectors
			@test typeof(det) == CylDetector
			@test typeof(det) != Detector
		end
		@test eltype(detectors) != CylDetector
		@test eltype(detectors) == Detector
		det1, det2, det3, det4 = detectors
		@test det1 <= det2 <= det3 <= det4

		detector_info_array = [5 1; 10 1; 15 1; 20 1//1]
		detectors = getDetectors(detector_info_array)
		for det = detectors
			@test typeof(det) == CylDetector
			@test typeof(det) != Detector
		end
		@test eltype(detectors) != CylDetector
		@test eltype(detectors) == Detector
		det1, det2, det3, det4 = detectors
		@test det1 <= det2 <= det3 <= det4

		detector_info_array = [5 1; 10 1; 15 1; 20 1.0]
		detectors = getDetectors(detector_info_array)
		for det = detectors
			@test typeof(det) == CylDetector
			@test typeof(det) != Detector
		end
		@test eltype(detectors) != CylDetector
		@test eltype(detectors) == Detector
		det1, det2, det3, det4 = detectors
		@test det1 <= det2 <= det3 <= det4

		detector_info_array = ["5" "0" "0" "0"; "10" "0" "0" "0"]
		@test_throws  MethodError getDetectors(detector_info_array; console_FB=false)
		detector_info_array = [5+1im 0 0 0; 5 10 0 0; 5 10 2 0; 5 10 2 5]
		@test_throws MethodError getDetectors(detector_info_array; console_FB=false)
		detector_info_array = detector_info_array = Matrix{Int}(undef, 0, 0)
		@test_throws ErrorException getDetectors(detector_info_array; console_FB=false)

		if true #!isapple()
			for q=["q", "Q"] , n=["", "n", "N", "fgdfgf", "qQ", "Qq"]
				write(stdin.buffer,"5\n" * "1\n" * "0\n" * "$q\n")
				@test getDetectors() == [Detector(5,1)]
				write(stdin.buffer,"5\n" * "1\n" * "0\n" * "$q\n")
				@test getDetectors(Matrix{Float64}(undef, 0, 0)) == [Detector(5,1)]

		 		write(stdin.buffer,"5\n" * "1\n" * "0\n" * "$n\n" *"55\n" * "11\n" * "0\n" * "$q\n")
				@test getDetectors(Matrix{Float64}(undef, 0, 0)) == [Detector(5, 1), Detector(55, 11)]
			end # for
		else
			@test_throws	ErrorException	write(stdin.buffer,"5\n" * "1\n" * "0\n" * "$q\n")
		end #if
    end # testset
println()
end # testset