text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Introduction to scikit-learn
#
# ## Basic preprocessing and model fitting
#
# In this notebook, we present how to build predictive models on tabular
# datasets.
#
# In particular we will highlight:
# * the difference between numerical and categorical variables;
# * the importance of scaling numerical variables;
# * typical ways to deal categorical variables;
# * how to train predictive models on different types of data;
# * evaluate the performance of a model via cross-validation.
#
# ## Introducing the dataset
#
# To this aim, we will use data from the 1994 Census bureau database. The goal
# with this data is to regress wages from heterogeneous data such as age,
# employment, education, family information, etc.
#
# Let's first load the data located in the `datasets` folder.

# %%
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")

# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

# %% [markdown]
# Let's have a look at the first records of this data frame:

# %%
df.head()

# %% [markdown]
# The target variable in our study will be the "class" column while we will use
# the other columns as input variables for our model. This target column divides
# the samples (also known as records) into two groups: high income (>50K) vs low
# income (<=50K). The resulting prediction problem is therefore a binary
# classification problem.
#
# For simplicity, we will ignore the "fnlwgt" (final weight) column that was
# crafted by the creators of the dataset when sampling the dataset to be
# representative of the full census database.

# %%
target_name = "class"
target = df[target_name].to_numpy()
target

# %%
data = df.drop(columns=[target_name, "fnlwgt"])
data.head()

# %% [markdown]
# We can check the number of samples and the number of features available in
# the dataset:

# %%
print(
    f"The dataset contains {data.shape[0]} samples and {data.shape[1]} "
    "features")

# %% [markdown]
# ## Working with numerical data
#
# Numerical data is the most natural type of data used in machine learning
# and can (almost) directly be fed to predictive models. We can quickly have a
# look at such data by selecting the subset of numerical columns from the
# original data.
#
# We will use this subset of data to fit a linear classification model to
# predict the income class.

# %%
data.columns

# %%
data.dtypes

# %%

# "i" denotes integer type, "f" denotes float type
numerical_columns = [
    c for c in data.columns if data[c].dtype.kind in ["i", "f"]]
numerical_columns

# %%
data_numeric = data[numerical_columns]
data_numeric.head()

# %% [markdown]
# When building a machine learning model, it is important to leave out a
# subset of the data which we can use later to evaluate the trained model.
# The data used to fit a model a called training data while the one used to
# assess a model are called testing data.
#
# Scikit-learn provides an helper function `train_test_split` which will
# split the dataset into a training and a testing set. It will ensure that
# the data are shuffled randomly before splitting the data.

# %%
from sklearn.model_selection import train_test_split

data_train, data_test, target_train, target_test = train_test_split(
    data_numeric, target, random_state=42)

print(
    f"The training dataset contains {data_train.shape[0]} samples and "
    f"{data_train.shape[1]} features")
print(
    f"The testing dataset contains {data_test.shape[0]} samples and "
    f"{data_test.shape[1]} features")

# %% [markdown]
# We will build a linear classification model called "Logistic Regression". The
# `fit` method is called to train the model from the input (features) and
# target data. Only the training data should be given for this purpose.
#
# In addition, check the time required to train the model and the number of
# iterations done by the solver to find a solution.
# %%
from sklearn.linear_model import LogisticRegression
import time

model = LogisticRegression(solver='lbfgs')
start = time.time()
model.fit(data_train, target_train)
elapsed_time = time.time() - start

print(f"The model {model.__class__.__name__} was trained in "
      f"{elapsed_time:.3f} seconds for {model.n_iter_} iterations")

# %% [markdown]
# Let's ignore the convergence warning for now and instead let's try
# to use our model to make some predictions on the first three records
# of the held out test set:

# %%
target_predicted = model.predict(data_test)
target_predicted[:5]

# %%
target_test[:5]

# %%
predictions = data_test.copy()
predictions['predicted-class'] = target_predicted
predictions['expected-class'] = target_test
predictions['correct'] = target_predicted == target_test
predictions.head()

# %% [markdown]
# To quantitatively evaluate our model, we can use the method `score`. It will
# compute the classification accuracy when dealing with a classification
# problem.

# %%
print(f"The test accuracy using a {model.__class__.__name__} is "
      f"{model.score(data_test, target_test):.3f}")

# %% [markdown]
# This is mathematically equivalent as computing the average number of time
# the model makes a correct prediction on the test set:

# %%
(target_test == target_predicted).mean()

# %% [markdown]
# ## Exercise 1
#
# - What would be the score of a model that always predicts `' >50K'`?
# - What would be the score of a model that always predicts `' <= 50K'`?
# - Is 81% or 82% accuracy a good score for this problem?
#
# Hint: You can compute the cross-validated of a
# [DummyClassifier](https://scikit-learn.org/stable/modules/model_evaluation.html#dummy-estimators)
# the performance of such baselines.
#
# Use the dedicated notebook to do this exercise.

# %% [markdown]
# Let's now consider the `ConvergenceWarning` message that was raised previously
# when calling the `fit` method to train our model. This warning informs us that
# our model stopped learning because it reached the maximum number of
# iterations allowed by the user. This could potentially be detrimental for the
# model accuracy. We can follow the (bad) advice given in the warning message
# and increase the maximum number of iterations allowed.

# %%
model = LogisticRegression(solver='lbfgs', max_iter=50000)
start = time.time()
model.fit(data_train, target_train)
elapsed_time = time.time() - start

# %%
print(
    f"The accuracy using a {model.__class__.__name__} is "
    f"{model.score(data_test, target_test):.3f} with a fitting time of "
    f"{elapsed_time:.3f} seconds in {model.n_iter_} iterations")

# %% [markdown]
# We now observe a longer training time but not significant improvement in
# the predictive performance. Instead of increasing the number of iterations, we
# can try to help fit the model faster by scaling the data first. A range of
# preprocessing algorithms in scikit-learn allows us to transform the input data
# before training a model. We can easily combine these sequential operations
# with a scikit-learn `Pipeline`, which chain together operations and can be
# used like any other classifier or regressor. The helper function
# `make_pipeline` will create a `Pipeline` by giving as arguments the successive
# transformations to perform followed by the classifier or regressor model.
#
# In our case, we will standardize the data and then train a new logistic
# regression model on that new version of the dataset set.

# %%
data_train.describe()

# %%
from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
data_train_scaled = scaler.fit_transform(data_train)
data_train_scaled

# %%
data_train_scaled = pd.DataFrame(data_train_scaled,
                                 columns=data_train.columns)
data_train_scaled.describe()

# %%
from sklearn.pipeline import make_pipeline

model = make_pipeline(StandardScaler(),
                      LogisticRegression(solver='lbfgs'))
start = time.time()
model.fit(data_train, target_train)
elapsed_time = time.time() - start

# %%
print(
    f"The accuracy using a {model.__class__.__name__} is "
    f"{model.score(data_test, target_test):.3f} with a fitting time of "
    f"{elapsed_time:.3f} seconds in {model[-1].n_iter_} iterations")

# %% [markdown]
# We can see that the training time and the number of iterations is much shorter
# while the predictive performance (accuracy) stays the same.
#
# In the previous example, we split the original data into a training set and a
# testing set. This strategy has several issues: in the setting where the amount
# of data is limited, the subset of data used to train or test will be small;
# and the splitting was done in a random manner and we have no information
# regarding the confidence of the results obtained.
#
# Instead, we can use cross-validation. Cross-validation consists of
# repeating this random splitting into training and testing sets and aggregating
# the model performance. By repeating the experiment, one can get an estimate of
# the variability of the model performance.
#
# The function `cross_val_score` allows for such experimental protocol by giving
# the model, the data and the target. Since there exists several
# cross-validation strategies, `cross_val_score` takes a parameter `cv` which
# defines the splitting strategy.

# %%
from sklearn.model_selection import cross_val_score

scores = cross_val_score(model, data_numeric, target, cv=5)
print(f"The different scores obtained are: \n{scores}")

# %%
print(f"The mean cross-validation accuracy is: "
      f"{scores.mean():.3f} +/- {scores.std():.3f}")

# %% [markdown]
# Note that by computing the standard-deviation of the cross-validation scores
# we can get an idea of the uncertainty of our estimation of the predictive
# performance of the model: in the above results, only the first 2 decimals seem
# to be trustworthy. Using a single train / test split would not allow us to
# know anything about the level of uncertainty of the accuracy of the model.
#
# Setting `cv=5` created 5 distinct splits to get 5 variations for the training
# and testing sets. Each training set is used to fit one model which is then
# scored on the matching test set. This strategy is called K-fold
# cross-validation where `K` corresponds to the number of splits.
#
# The following matplotlib code helps visualize how the dataset is partitioned
# into train and test samples at each iteration of the cross-validation
# procedure:

# %%
# %matplotlib inline
from sklearn.model_selection import KFold
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Patch

cmap_cv = plt.cm.coolwarm

def plot_cv_indices(cv, X, y, ax, lw=20):
    """Create a sample plot for indices of a cross-validation object."""
    splits = list(cv.split(X=X, y=y))
    n_splits = len(splits)

    # Generate the training/testing visualizations for each CV split
    for ii, (train, test) in enumerate(splits):
        # Fill in indices with the training/test groups
        indices = np.zeros(shape=X.shape[0], dtype=np.int32)
        indices[train] = 1

        # Visualize the results
        ax.scatter(range(len(indices)), [ii + .5] * len(indices),
                   c=indices, marker='_', lw=lw, cmap=cmap_cv,
                   vmin=-.2, vmax=1.2)

    # Formatting
    yticklabels = list(range(n_splits))
    ax.set(yticks=np.arange(n_splits + 2) + .5,
           yticklabels=yticklabels, xlabel='Sample index',
           ylabel="CV iteration", ylim=[n_splits + .2,
                                        -.2], xlim=[0, 100])
    ax.set_title('{}'.format(type(cv).__name__), fontsize=15)
    return ax

# %%
# Some random data points
n_points = 100
X = np.random.randn(n_points, 10)
y = np.random.randn(n_points)

fig, ax = plt.subplots(figsize=(10, 6))
cv = KFold(5)
_ = plot_cv_indices(cv, X, y, ax)

# TODO: add summary here


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.3.0
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Introduction to scikit-learn: basic model hyper-parameters tuning
#
# The process of learning a predictive model is driven by a set of internal
# parameters and a set of training data. These internal parameters are called
# hyper-parameters and are specific for each family of models. In addition, a
# specific set of parameters are optimal for a specific dataset and thus they
# need to be optimized.
#
# This notebook shows:
# * the influence of changing model parameters;
# * how to tune these hyper-parameters;
# * how to evaluate the model performance together with hyper-parameter
#   tuning.

# %%
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")
# Or use the local copy:
# df = pd.read_csv(os.path.join("..", "datasets", "adult-census.csv"))

# %%
target_name = "class"
target = df[target_name].to_numpy()
target

# %%
data = df.drop(columns=[target_name, "fnlwgt"])
data.head()

# %% [markdown]
# Once the dataset is loaded, we split it into a training and testing sets.

# %%
from sklearn.model_selection import train_test_split

df_train, df_test, target_train, target_test = train_test_split(
    data, target, random_state=42)

# %% [markdown]
# Then, we define the preprocessing pipeline to transform differently
# the numerical and categorical data.

# %%
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OrdinalEncoder

categorical_columns = [
    'workclass', 'education', 'marital-status', 'occupation',
    'relationship', 'race', 'native-country', 'sex']

categories = [
    data[column].unique() for column in data[categorical_columns]]

categorical_preprocessor = OrdinalEncoder(categories=categories)

preprocessor = ColumnTransformer([
    ('cat-preprocessor', categorical_preprocessor,
     categorical_columns),], remainder='passthrough',
                                 sparse_threshold=0)

# %% [markdown]
# Finally, we use a tree-based classifier (i.e. histogram gradient-boosting) to
# predict whether or not a person earns more than 50,000 dollars a year.

# %%
# %%time
# for the moment this line is required to import HistGradientBoostingClassifier
from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.pipeline import Pipeline

model = Pipeline([
    ("preprocessor", preprocessor),
    ("classifier",
     HistGradientBoostingClassifier(max_leaf_nodes=16,
                                    learning_rate=0.05,
                                    random_state=42)),])
model.fit(df_train, target_train)

print(
    f"The test accuracy score of the gradient boosting pipeline is: "
    f"{model.score(df_test, target_test):.2f}")

# %% [markdown]
# ## Quizz
#
# 1. What is the default value of the `learning_rate` parameter of the `HistGradientBoostingClassifier` class? ([link to the API documentation](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.HistGradientBoostingClassifier.html#sklearn-ensemble-histgradientboostingclassifier))
#
# 2. Try to edit the code of the previous cell to set the learning rate parameter to 10. Does this increase the accuracy of the model?
#
# 3. Decrease progressively value of `learning_rate`: can you find a value that yields an accuracy higher than 0.87?
#
# 4. Fix `learning_rate` to 0.05 and try setting the value of `max_leaf_nodes` to the minimum value of 2. Does not improve the accuracy?
#
# 5. Try to progressively increase the value of `max_leaf_nodes` to 256 by taking powers of 2. What do you observe?

# %% [markdown]
# ## The issue of finding the best model parameters
#
# In the previous example, we created an histogram gradient-boosting classifier
# using the default parameters by omitting to explicitely set these parameters.
#
# However, there is no reasons that this set of parameters are optimal for our
# dataset. For instance, fine-tuning the histogram gradient-boosting can be
# achieved by finding the best combination of the following parameters: (i)
# `learning_rate`, (ii) `min_samples_leaf`, and (iii) `max_leaf_nodes`.
# Nevertheless, finding this combination manually will be tedious. Indeed,
# there are relationship between these parameters which are difficult to find
# manually: increasing the depth of trees (increasing `max_samples_leaf`)
# should be associated with a lower learning-rate.
#
# Scikit-learn provides tools to explore and evaluate the parameters
# space.
# %% [markdown]
# ## Finding the best model hyper-parameters via exhaustive parameters search
#
# Our goal is to find the best combination of the parameters stated above.
#
# In short, we will set these parameters with some defined values, train our
# model on some data, and evaluate the model performance on some left out data.
# Ideally, we will select the parameters leading to the optimal performance on
# the testing set.

# %% [markdown]
# The first step is to find the name of the parameters to be set. We use the
# method `get_params()` to get this information. For instance, for a single
# model like the `HistGradientBoostingClassifier`, we can get the list such as:

# %%
print("The hyper-parameters are for a histogram GBDT model are:")
for param_name in HistGradientBoostingClassifier().get_params().keys(
):
    print(param_name)

# %% [markdown]
# When the model of interest is a `Pipeline`, i.e. a serie of transformers and
# a predictor, the name of the estimator will be added at the front of the
# parameter name with a double underscore ("dunder") in-between (e.g.
# `estimator__parameters`).

# %%
print("The hyper-parameters are for the full-pipeline are:")
for param_name in model.get_params().keys():
    print(param_name)

# %% [markdown]
# The parameters that we want to set are:
# - `'classifier__learning_rate'`: this parameter will
#   control the ability of a new tree to correct the error of the previous
#   sequence of trees;
# - `'classifier__max_leaf_nodes'`: this parameter will
#   control the depth of each tree.

# %% [markdown]
# ## Exercises:
#
# Use the previously defined model (called `model`) and using two nested `for`
# loops, make a search of the best combinations of the `learning_rate` and
# `max_leaf_nodes` parameters. In this regard, you will need to train and test
# the model by setting the parameters. The evaluation of the model should be
# performed using `cross_val_score`. We can propose to define the following
# parameters search:
# - `learning_rate` for the values 0.01, 0.1, and 1;
# - `max_leaf_nodes` for the values 5, 25, 45.

# %% [markdown]
# ## Automated parameter tuning via grid-search
#
# Instead of manually writting the two `for` loops, scikit-learn provides a
# class called `GridSearchCV` which implement the exhaustive search implemented
# during the exercise.
#
# Let see how to use the `GridSearchCV` estimator for doing such search.
# Since the grid-search will be costly, we will only explore the combination
# learning-rate and the maximum number of nodes.

# %%
# %%time
import numpy as np
from sklearn.model_selection import GridSearchCV

param_grid = {
    'classifier__learning_rate': (0.05, 0.1, 0.5, 1, 5),
    'classifier__max_leaf_nodes': (3, 10, 30, 100),}
model_grid_search = GridSearchCV(model, param_grid=param_grid,
                                 n_jobs=4, cv=2)
model_grid_search.fit(df_train, target_train)

print(f"The test accuracy score of the grid-searched pipeline is: "
      f"{model_grid_search.score(df_test, target_test):.2f}")

# %% [markdown]
# The `GridSearchCV` estimator takes a `param_grid` parameter which defines
# all hyper-parameters and their associated values. The grid-search will be in
# charge of creating all possible combinations and test them.
#
# The number of combinations will be equal to the cardesian product of the
# number of values to explore for each parameter (e.g. in our example 3 x 3
# combinations). Thus, adding new parameters with their associated values to be
# explored become rapidly computationally expensive.
#
# Once the grid-search is fitted, it can be used as any other predictor by
# calling `predict` and `predict_proba`. Internally, it will use the model with
# the best parameters found during `fit`.
#
# Get predictions for the 5 first samples using the estimator with the best
# parameters.

# %%
model_grid_search.predict(df_test.iloc[0:5])

# %% [markdown]
# You can know about these parameters by looking at the `best_params_`
# attribute.

# %%
print(f"The best set of parameters is: "
      f"{model_grid_search.best_params_}")

# %% [markdown]
# In addition, we can inspect all results which are stored in the attribute
# `cv_results_` of the grid-search. We will filter some specific columns to
# from these results

# %%
cv_results = pd.DataFrame(model_grid_search.cv_results_).sort_values(
    "mean_test_score", ascending=False)
cv_results.head()

# %% [markdown]
# Let us focus on the most interesting columns and shorten the parameter names to remove the `"param_classifier__"` prefix for readability:

# %%
# get the parameter names
column_results = [f"param_{name}" for name in param_grid.keys()]
column_results += [
    "mean_test_score", "std_test_score", "rank_test_score"]
cv_results = cv_results[column_results]

# %%
def shorten_param(param_name):
    if "__" in param_name:
        return param_name.rsplit("__", 1)[1]
    return param_name

cv_results = cv_results.rename(shorten_param, axis=1)
cv_results

# %% [markdown]
# With only 2 parameters, we might want to visualize the grid-search as a
# heatmap. We need to transform our `cv_results` into a dataframe where the
# rows will correspond to the learning-rate values and the columns will
# correspond to the maximum number of leaf and the content of the dataframe
# will be the mean test scores.

# %%
pivoted_cv_results = cv_results.pivot_table(
    values="mean_test_score", index=["learning_rate"],
    columns=["max_leaf_nodes"])

pivoted_cv_results

# %%
import matplotlib.pyplot as plt
from seaborn import heatmap

ax = heatmap(pivoted_cv_results, annot=True, cmap="YlGnBu", vmin=0.7,
             vmax=0.9)
ax.invert_yaxis()

# %% [markdown]
# The above tables highlights the following things:
#
# - for too high values of the value of `learning_rate`, the performance of the model is degraded and adjusting the value of `max_leaf_nodes` cannot fix that problem;
# - outside of this pathological region, we observe that the optimal choice of `max_leaf_nodes` depends on the value of `learning_rate`;
# - in particular, we observe a "diagonal" of good models with an accuracy close to the maximal of 0.87: when the value of `max_leaf_nodes` is increased, one should increase the value of `learning_rate` accordingly to preserve a good accuracy.
#
# The precise meaning of those two parameters will be explained in a latter notebook.
#
# For now we will note that, in general, **there is no unique optimal parameter setting**: 6 models out of the 16 parameter configuration reach the maximal accuracy (up to smal random fluctuations caused by the sampling of the training set).

# %% [markdown]
# ## Hyper-parameter tuning with Random Search
#
#
# With the `GridSearchCV` estimator, the parameters need to be specified
# explicitely. We mentioned that exploring a large number of values for
# different parameters will be quickly untractable.
#
# Instead, we can randomly generate the parameter candidates. The
# `RandomSearchCV` allows for such stochastic search. It is used similarly to
# the `GridSearchCV` but the sampling distributions need to be specified
# instead of the parameter values. For instance, we will draw candidates using
# a log-uniform distribution also called reciprocal distribution. In addition,
# we will optimize 2 other parameters:
# - `max_iter`: it corresponds to the number of trees in the ensemble;
# - `min_samples_leaf`: it corresponds to the minimum number of samples
#   required in a leaf.

# %%
from scipy.stats import reciprocal
from sklearn.model_selection import RandomizedSearchCV
from pprint import pprint

class reciprocal_int:
    """Integer valued version of the log-uniform distribution"""
    def __init__(self, a, b):
        self._distribution = reciprocal(a, b)

    def rvs(self, *args, **kwargs):
        """Random variable sample"""
        return self._distribution.rvs(*args, **kwargs).astype(int)

param_distributions = {
    'classifier__l2_regularization': reciprocal(1e-6, 1e3),
    'classifier__learning_rate': reciprocal(0.001, 10),
    'classifier__max_leaf_nodes': reciprocal_int(2, 256),
    'classifier__min_samples_leaf': reciprocal_int(1, 100),
    'classifier__max_bins': reciprocal_int(2, 255),}
model_random_search = RandomizedSearchCV(
    model, param_distributions=param_distributions, n_iter=10,
    n_jobs=4, cv=5)
model_random_search.fit(df_train, target_train)

print(f"The test accuracy score of the best model is "
      f"{model_random_search.score(df_test, target_test):.2f}")

# %%
print("The best parameters are:")
pprint(model_random_search.best_params_)

# %% [markdown]
# We can inspect the results using the attributes `cv_results` as we previously
# did.

# %%
# get the parameter names
column_results = [
    f"param_{name}" for name in param_distributions.keys()]
column_results += [
    "mean_test_score", "std_test_score", "rank_test_score"]

cv_results = pd.DataFrame(model_random_search.cv_results_)
cv_results = cv_results[column_results].sort_values(
    "mean_test_score", ascending=False)
cv_results = cv_results.rename(shorten_param, axis=1)
cv_results

# %% [markdown]
# In practice, a randomized hyper-parameter search is usually run with a large number of
# iterations. In order to avoid the computation cost and still make a decent
# analysis, we load the results obtained from a similar search with 200
# iterations.

# %%
# model_random_search = RandomizedSearchCV(
#     model, param_distributions=param_distributions, n_iter=500,
#     n_jobs=4, cv=5)
# model_random_search.fit(df_train, target_train)
# cv_results =  pd.DataFrame(model_random_search.cv_results_)
# cv_results.to_csv("../figures/randomized_search_results.csv")

# %%
cv_results = pd.read_csv("../figures/randomized_search_results.csv",
                         index_col=0)

# %% [markdown]
# As we have more than 2 paramters in our grid-search, we cannot visualize the
# results using a heatmap. However, we can us a parallel coordinates plot.

# %%
(cv_results[column_results].rename(
    shorten_param, axis=1).sort_values("mean_test_score"))

# %%
import plotly.express as px

fig = px.parallel_coordinates(
    cv_results.rename(shorten_param, axis=1).apply({
        "learning_rate": np.log10,
        "max_leaf_nodes": np.log2,
        "max_bins": np.log2,
        "min_samples_leaf": np.log10,
        "l2_regularization": np.log10,
        "mean_test_score": lambda x: x,}),
    color="mean_test_score",
    color_continuous_scale=px.colors.sequential.Viridis,
)
fig.show()

# %% [markdown]
# The parallel coordinates plot will display the values of the hyper-parameters
# on different columns while the performance metric is color coded. Thus, we
# are able to quickly inspect if there is a range of hyper-parameters which is
# working or not.
#
# Note that we **transformed most axis values by taking a log10 or log2** to
# spread the active ranges and improve the readability of the plot.
#
# It is possible to **select a range of results by clicking and holding on
# any axis** of the parallel coordinate plot. You can then slide (move)
# the range selection and cross two selections to see the intersections.

# %% [markdown]
# **Quizz**
#
#
# Select the worst performing models (for instance models with a "mean_test_score" lower than 0.7): what do have all these moels in common (choose one):
#
#
# |                               |      |
# |-------------------------------|------|
# | too large `l2_regularization` |      |
# | too small `l2_regularization` |      |
# | too large `learning_rate`     |      |
# | too low `learning_rate`       |      |
# | too large `max_bins`          |      |
# | too large `max_bins`          |      |
#
#
# Using the above plot, identify ranges of values for hyperparameter that always prevent the model to reach a test score higher than 0.86, irrespective of the other values:
#
#
# |                               | True | False |
# |-------------------------------|------|-------|
# | too large `l2_regularization` |      |       |
# | too small `l2_regularization` |      |       |
# | too large `learning_rate`     |      |       |
# | too low `learning_rate`       |      |       |
# | too large `max_bins`          |      |       |
# | too large `max_bins`          |      |       |

# %% [markdown]
# ## Exercises:
#
# - Build a machine learning pipeline:
#   * preprocess the categorical columns using a `OneHotEncoder` and use
#     a `StandardScaler` to normalize the numerical data.
#   * use a `LogisticRegression` as a predictive model.
# - Make an hyper-parameters search using `RandomizedSearchCV` and tuning the
#   parameters:
#   * `C` with values ranging from 0.001 to 10. You can use a reciprocal
#     distribution (i.e. `scipy.stats.reciprocal`);
#   * `solver` with possible values being `"liblinear"` and `"lbfgs"`;
#   * `penalty` with possible values being `"l2"` and `"l1"`;
#   * `drop` with possible values being `None` or `"first"`.
#
# You might get some `FitFailedWarning` and try to explain why.

# %% [markdown]
# ## Combining evaluation and hyper-parameters search
#
# Cross-validation was used for searching for the best model parameters. We
# previously evaluated model performance through cross-validation as well. If
# we would like to combine both aspects, we need to perform a **"nested"
# cross-validation**. The "outer" cross-validation is applied to assess the model
# while the "inner" cross-validation sets the hyper-parameters of the model on
# the data set provided by the "outer" cross-validation.
#
#
# In practice, it can be implemented by calling `cross_val_score` or
# `cross_validate` on an instance of `GridSearchCV`, `RandomSearchCV`, or any
# other `EstimatorCV` class.

# %%
from sklearn.model_selection import cross_val_score

# recall the definition of our grid-search
param_distributions = {
    'classifier__max_iter': reciprocal_int(10, 50),
    'classifier__learning_rate': reciprocal(0.01, 10),
    'classifier__max_leaf_nodes': reciprocal_int(2, 16),
    'classifier__min_samples_leaf': reciprocal_int(1, 50),}
model_random_search = RandomizedSearchCV(
    model, param_distributions=param_distributions, n_iter=10,
    n_jobs=4, cv=5)

scores = cross_val_score(model_random_search, data, target, n_jobs=4,
                         cv=5)

# %%
print(f"The cross-validated accuracy score is:"
      f" {scores.mean():.3f} +- {scores.std():.3f}")

# %%
print("The scores obtained for each CV split are:")
print(scores)

# %% [markdown]
# Be aware that the best model found for each split of the outer cross-validation loop might not share the same hyper-parameter values.
#
# When analyzing such model, you should not only look at the
# overall model performance but look at the hyper-parameters variations as
# well.

# %% [markdown]
# ## In this notebook, we have:
#
# * manually tuned the hyper-parameters of a machine-learning pipeline;
# * automatically tuned the hyper-parameters of a machine-learning pipeline by
#   by exhaustively searching the best combination of parameters from a defined
#   grid;
# * automatically tuned the hyper-parameters of a machine-learning pipeline by
#   drawing values candidates from some predefined distributions;
# * nested an hyper-parameters tuning procedure within an cross-validation
#   evaluation procedure.
#
# ## Main take-away points
#
# * a grid-search is a costly exhaustive search and does scale with the number of
#   parameters to search;
# * a randomized-search will always run with a fixed given budget;
# * when assessing the performance of a model, hyper-parameters search should
#   be tuned on the training data of a predifined train test split;
# * alternatively it is possible to nest parameter tuning within a
#   cross-validation scheme.

# %%

repo_name: lesteve/scikit-learn-tutorial

repo_language: Jupyter Notebook


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Exercise 03
#
# The goal of this exercise is to evaluate the impact of feature preprocessing on a pipeline that uses a  decision-tree-based classifier instead of logistic regression.
#
# - The first question is to empirically evaluate whether scaling numerical feature is helpful or not;
#
# - The second question is to evaluate whether it is empirically better (both from a computational and a statistical perspective) to use integer coded or one-hot encoded categories.

# %%
import pandas as pd
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OrdinalEncoder
from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")

# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

# %%
target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=[target_name, "fnlwgt"])

# %%
numerical_columns = [
    c for c in data.columns if data[c].dtype.kind in ["i", "f"]]
categorical_columns = [
    c for c in data.columns if data[c].dtype.kind not in ["i", "f"]]

categories = [
    data[column].unique() for column in data[categorical_columns]]

# %% [markdown]
# ## Reference pipeline (no numerical scaling and integer-coded categories)
#
# First let's time the pipeline we used in the main notebook to serve as a reference:

# %%
# %%time

preprocessor = ColumnTransformer([
    ('categorical', OrdinalEncoder(categories=categories),
     categorical_columns),], remainder="passthrough")

model = make_pipeline(preprocessor, HistGradientBoostingClassifier())
scores = cross_val_score(model, data, target)
print(f"The different scores obtained are: \n{scores}")
print(f"The accuracy is: {scores.mean():.3f} +- {scores.std():.3f}")

# %% [markdown]
# ## Scaling numerical features
#
# Let's write a similar pipeline that also scales the numerical features using `StandardScaler` (or similar):

# %%
# TODO write me!

# %% [markdown]
# ## One-hot encoding of categorical variables
#
# For linear models, we have observed that integer coding of categorical
# variables can be very detrimental. However for
# `HistGradientBoostingClassifier` models, it does not seem to be the
# case as the cross-validation of the reference pipeline with
# `OrdinalEncoder` is good.
#
# Let's see if we can get an even better accuracy with `OneHotEncoding`.
#
# Hint: `HistGradientBoostingClassifier` does not yet support sparse input data. You might want to use
# `OneHotEncoder(handle_unknown="ignore", sparse=False)` to force the use a dense representation as a workaround.

# %%
# TODO: write me!


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# #  Exercise 01
#
# The goal of is to compare the performance of our classifier (81% accuracy) to some baseline classifiers that  would ignore the input data and instead make constant predictions.
#
# The online [documentation for DummyClassifier](https://scikit-learn.org/stable/modules/model_evaluation.html#dummy-estimators) gives instructions on how to use it.

# %%
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")

# %%
target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=[target_name, "fnlwgt"])
numerical_columns = [
    c for c in data.columns if data[c].dtype.kind in ["i", "f"]]
data_numeric = data[numerical_columns]

# %%
from sklearn.model_selection import cross_val_score
from sklearn.dummy import DummyClassifier

# TODO: write me!


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Exercise 02
#
# The goal of this exercise is to evalutate the impact of using an arbitrary
# integer encoding for categorical variables along with a linear
# classification model such as Logistic Regression.
#
# To do so, let's try to use `OrdinalEncoder` to preprocess the categorical
# variables. This preprocessor is assembled in a pipeline with
# `LogisticRegression`. The performance of the pipeline can be evaluated as
# usual by cross-validation and then compared to the score obtained when using
# `OneHotEncoding` or to some other baseline score.
#
# Because `OrdinalEncoder` can raise errors if it sees an unknown category at
# prediction time, we need to pre-compute the list of all possible categories
# ahead of time:
#
# ```python
# categories = [data[column].unique()
#               for column in data[categorical_columns]]
# OrdinalEncoder(categories=categories)
# ```

# %%
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")

# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

# %%
target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=[target_name, "fnlwgt"])
categorical_columns = [
    c for c in data.columns if data[c].dtype.kind not in ["i", "f"]]
data_categorical = data[categorical_columns]

# %%
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OrdinalEncoder
from sklearn.linear_model import LogisticRegression

# TODO: write me!


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# #  Exercise 01
# The goal is to write an exhaustive search to find the best parameters
# combination maximizing the model performance

# %%
import pandas as pd

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import OrdinalEncoder
from sklearn.model_selection import RandomizedSearchCV
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
# This line is currently required to import HistGradientBoostingClassifier
from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier

from scipy.stats import expon, uniform
from scipy.stats import randint

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")
# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=target_name)

df_train, df_test, target_train, target_test = train_test_split(
    data, target, random_state=42)

from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OrdinalEncoder

categorical_columns = [
    'workclass', 'education', 'marital-status', 'occupation',
    'relationship', 'race', 'native-country', 'sex']

categories = [data[column].unique()
              for column in data[categorical_columns]]

categorical_preprocessor = OrdinalEncoder(categories=categories)

preprocessor = ColumnTransformer(
    [('cat-preprocessor', categorical_preprocessor, categorical_columns)],
    remainder='passthrough', sparse_threshold=0)

from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.pipeline import make_pipeline

model = make_pipeline(
    preprocessor, HistGradientBoostingClassifier(random_state=42))

# %% [markdown]
# TODO: write your solution here
#
# Use the previously defined model (called `model`) and using two nested `for`
# loops, make a search of the best combinations of the `learning_rate` and
# `max_leaf_nodes` parameters. In this regard, you will need to train and test
# the model by setting the parameters. The evaluation of the model should be
# performed using `cross_val_score`. We can propose to define the following
# parameters search:
# - `learning_rate` for the values 0.01, 0.1, and 1;
# - `max_leaf_nodes` for the values 5, 25, 45.


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Solution for Exercise 02
#
# The goal of this exercise is to evalutate the impact of using an arbitrary
# integer encoding for categorical variables along with a linear
# classification model such as Logistic Regression.
#
# To do so, let's try to use `OrdinalEncoder` to preprocess the categorical
# variables. This preprocessor is assembled in a pipeline with
# `LogisticRegression`. The performance of the pipeline can be evaluated as
# usual by cross-validation and then compared to the score obtained when using
# `OneHotEncoding` or to some other baseline score.
#
# Because `OrdinalEncoder` can raise errors if it sees an unknown category at
# prediction time, we need to pre-compute the list of all possible categories
# ahead of time:
#
# ```python
# categories = [data[column].unique()
#               for column in data[categorical_columns]]
# OrdinalEncoder(categories=categories)
# ```

# %%
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")

# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

# %%
target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=[target_name, "fnlwgt"])
categorical_columns = [
    c for c in data.columns if data[c].dtype.kind not in ["i", "f"]]
data_categorical = data[categorical_columns]

# %%
categories = [
    data[column].unique() for column in data[categorical_columns]]

categories

# %%
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OrdinalEncoder
from sklearn.linear_model import LogisticRegression

model = make_pipeline(
    OrdinalEncoder(categories=categories),
    LogisticRegression(solver='lbfgs', max_iter=1000))
scores = cross_val_score(model, data_categorical, target)
print(f"The different scores obtained are: \n{scores}")

# %%
print(f"The accuracy is: {scores.mean():.3f} +- {scores.std():.3f}")

# %% [markdown]
# Using an arbitrary mapping from string labels to integers as done here causes the linear model to make bad assumptions on the relative ordering of  categories.
#
# This prevent the model to learning anything predictive enough and the cross-validated score is even lower that the baseline we obtained by ignoring the input data and just always predict the most frequent class:

# %%
from sklearn.dummy import DummyClassifier

scores = cross_val_score(DummyClassifier(strategy="most_frequent"),
                         data_categorical, target)
print(f"The different scores obtained are: \n{scores}")
print(f"The accuracy is: {scores.mean():.3f} +- {scores.std():.3f}")

# %% [markdown]
# By comparison, a categorical encoding that does not assume any ordering in the
# categories can lead to a significantly higher score:

# %%
from sklearn.preprocessing import OneHotEncoder

model = make_pipeline(
    OneHotEncoder(handle_unknown="ignore"),
    LogisticRegression(solver='lbfgs', max_iter=1000))
scores = cross_val_score(model, data_categorical, target)
print(f"The different scores obtained are: \n{scores}")
print(f"The accuracy is: {scores.mean():.3f} +- {scores.std():.3f}")


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Exercise 02
# The goal is to find the best set of hyper-parameters which maximize the
# performance on a training set.

# %%
import numpy as np
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")
# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=target_name)

from sklearn.model_selection import train_test_split

df_train, df_test, target_train, target_test = train_test_split(
    data, target, random_state=42)

# %% [markdown]
# TODO: create your machine learning pipeline
#
# You should:
# * preprocess the categorical columns using a `OneHotEncoder` and use a
#   `StandardScaler` to normalize the numerical data.
# * use a `LogisticRegression` as a predictive model.

# %% [markdown]
# Start by defining the columns and the preprocessing pipelines to be applied
# on each columns.
# %%

from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import StandardScaler

# %% [markdown]
# Subsequently, create a `ColumnTransformer` to redirect the specific columns
# a preprocessing pipeline.
# %%

from sklearn.compose import ColumnTransformer

# %% [markdown]
# Finally, concatenate the preprocessing pipeline with a logistic regression.
# %%

from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression

# %% [markdown]
# TODO: make your random search
#
# Use a `RandomizedSearchCV` to find the best set of hyper-parameters by tuning
# the following parameters for the `LogisticRegression` model:
# - `C` with values ranging from 0.001 to 10. You can use a reciprocal
#   distribution (i.e. `scipy.stats.reciprocal`);
# - `solver` with possible values being `"liblinear"` and `"lbfgs"`;
# - `penalty` with possible values being `"l2"` and `"l1"`;
# In addition, try several preprocessing strategies with the `OneHotEncoder`
# by always (or not) dropping the first column when encoding the categorical
# data.
#
# Notes: You can accept failure during a grid-search or a randomized-search
# by settgin `error_score` to `np.nan` for instance.


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Solution for Exercise 01
#
# The goal of is to compare the performance of our classifier to some baseline classifier that  would ignore the input data and instead make constant predictions:

# %%
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")

# %%
target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=[target_name, "fnlwgt"])
numerical_columns = [
    c for c in data.columns if data[c].dtype.kind in ["i", "f"]]
data_numeric = data[numerical_columns]

# %%
from sklearn.model_selection import cross_val_score
from sklearn.dummy import DummyClassifier

high_revenue_clf = DummyClassifier(strategy="constant",
                                   constant=" >50K")
scores = cross_val_score(high_revenue_clf, data_numeric, target)
print(f"{scores.mean():.3f} +/- {scores.std():.3f}")

# %%
low_revenue_clf = DummyClassifier(strategy="constant",
                                  constant=" <=50K")
scores = cross_val_score(low_revenue_clf, data_numeric, target)
print(f"{scores.mean():.3f} +/- {scores.std():.3f}")

# %%
most_freq_revenue_clf = DummyClassifier(strategy="most_frequent")
scores = cross_val_score(most_freq_revenue_clf, data_numeric, target)
print(f"{scores.mean():.3f} +/- {scores.std():.3f}")

# %% [markdown]
# So 81% accuracy is significantly better than 76% which is the score of a baseline model that would always predict the most frequent class which is the low revenue class: `" <=50K"`.
#
# In this dataset, we can see that the target classes are imbalanced: almost 3/4 of the records are people with a revenue below 50K:

# %%
df["class"].value_counts()

# %%
(target == " <=50K").mean()


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Solution for Exercise 03
#
# The goal of this exercise is to evaluate the impact of feature preprocessing on a pipeline that uses a  decision-tree-based classifier instead of logistic regression.
#
# - The first question is to empirically evaluate whether scaling numerical feature is helpful or not;
#
# - The second question is to evaluate whether it is empirically better (both from a computational and a statistical perspective) to use integer coded or one-hot encoded categories.
#
#
# Hint: `HistGradientBoostingClassifier` does not yet support sparse input data. You might want to use
# `OneHotEncoder(handle_unknown="ignore", sparse=False)` to force the use a dense representation as a workaround.

# %%
import pandas as pd
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OrdinalEncoder
from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")

# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

# %%
target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=[target_name, "fnlwgt"])

# %%
numerical_columns = [
    c for c in data.columns if data[c].dtype.kind in ["i", "f"]]
categorical_columns = [
    c for c in data.columns if data[c].dtype.kind not in ["i", "f"]]

categories = [
    data[column].unique() for column in data[categorical_columns]]

# %% [markdown]
# ## Reference pipeline (no numerical scaling and integer-coded categories)
#
# First let's time the pipeline we used in the main notebook to serve as a reference:

# %%
# %%time

preprocessor = ColumnTransformer([
    ('categorical', OrdinalEncoder(categories=categories),
     categorical_columns),], remainder="passthrough")

model = make_pipeline(preprocessor, HistGradientBoostingClassifier())
scores = cross_val_score(model, data, target)
print(f"The different scores obtained are: \n{scores}")
print(f"The accuracy is: {scores.mean():.3f} +- {scores.std():.3f}")

# %% [markdown]
# ## Scaling numerical features

# %%
# %%time
from sklearn.preprocessing import StandardScaler

preprocessor = ColumnTransformer([
    ('numerical', StandardScaler(), numerical_columns),
    ('categorical', OrdinalEncoder(categories=categories),
     categorical_columns),])

model = make_pipeline(preprocessor, HistGradientBoostingClassifier())
scores = cross_val_score(model, data, target)
print(f"The different scores obtained are: \n{scores}")
print(f"The accuracy is: {scores.mean():.3f} +- {scores.std():.3f}")

# %% [markdown]
# ### Analysis
#
# We can observe that both the accuracy and the training time are approximately the same as the reference pipeline (any time difference you might observe is not significant).
#
# Scaling numerical features is indeed useless for most decision tree models in general and for `HistGradientBoostingClassifier` in particular.

# %% [markdown]
# ## One-hot encoding of categorical variables
#
# For linear models, we have observed that integer coding of categorical
# variables can be very detrimental. However for
# `HistGradientBoostingClassifier` models, it does not seem to be the
# case as the cross-validation of the reference pipeline with
# `OrdinalEncoder` is good.
#
# Let's see if we can get an even better accuracy with `OneHotEncoding`:

# %%
# %%time
from sklearn.preprocessing import OneHotEncoder

preprocessor = ColumnTransformer([
    ('categorical',
     OneHotEncoder(handle_unknown="ignore", sparse=False),
     categorical_columns),], remainder="passthrough")

model = make_pipeline(preprocessor, HistGradientBoostingClassifier())
scores = cross_val_score(model, data, target)
print(f"The different scores obtained are: \n{scores}")
print(f"The accuracy is: {scores.mean():.3f} +- {scores.std():.3f}")

# %% [markdown]
# ### Analysis
#
# From an accuracy point of view, the result is almost exactly the same.
# The reason is that `HistGradientBoostingClassifier` is expressive
# and robust enough to deal with misleading ordering of integer coded
# categories (which was not the case for linear models).
#
# However from a computation point of view, the training time is
# significantly longer: this is caused by the fact that `OneHotEncoder`
# generates approximately 10 times more features than `OrdinalEncoder`.
#
# Note that the current implementation `HistGradientBoostingClassifier`
# is still incomplete, and once sparse representation are handled
# correctly, training time might improve with such kinds of encodings.
#
# The main take away message is that arbitrary integer coding of
# categories is perfectly fine for `HistGradientBoostingClassifier`
# and yields fast training times.


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Exercise 02
# The goal is to find the best set of hyper-parameters which maximize the
# performance on a training set.

# %%
import numpy as np
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")
# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=target_name)

from sklearn.model_selection import train_test_split

df_train, df_test, target_train, target_test = train_test_split(
    data, target, random_state=42)

# %% [markdown]
# TODO: create your machine learning pipeline
#
# You should:
# * preprocess the categorical columns using a `OneHotEncoder` and use a
#   `StandardScaler` to normalize the numerical data.
# * use a `LogisticRegression` as a predictive model.

# %% [markdown]
# Start by defining the columns and the preprocessing pipelines to be applied
# on each columns.
# %%
categorical_columns = [
    'workclass', 'education', 'marital-status', 'occupation',
    'relationship', 'race', 'native-country', 'sex']

categories = [data[column].unique()
              for column in data[categorical_columns]]

numerical_columns = [
    'age', 'capital-gain', 'capital-loss', 'hours-per-week']

from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import StandardScaler

categorical_processor = OneHotEncoder(categories=categories)
numerical_processor = StandardScaler()

# %% [markdown]
# Subsequently, create a `ColumnTransformer` to redirect the specific columns
# a preprocessing pipeline.
# %%

from sklearn.compose import ColumnTransformer

preprocessor = ColumnTransformer(
    [('cat-preprocessor', categorical_processor, categorical_columns),
     ('num-preprocessor', numerical_processor, numerical_columns)]
)

# %% [markdown]
# Finally, concatenate the preprocessing pipeline with a logistic regression.
# %%

from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression

model = make_pipeline(preprocessor, LogisticRegression())

# %% [markdown]
# Use a `RandomizedSearchCV` to find the best set of hyper-parameters by tuning
# the following parameters for the `LogisticRegression` model:
# - `C` with values ranging from 0.001 to 10. You can use a reciprocal
#   distribution (i.e. `scipy.stats.reciprocal`);
# - `solver` with possible values being `"liblinear"` and `"lbfgs"`;
# - `penalty` with possible values being `"l2"` and `"l1"`;
# In addition, try several preprocessing strategies with the `OneHotEncoder`
# by always (or not) dropping the first column when encoding the categorical
# data.
#
# Notes: You can accept failure during a grid-search or a randomized-search
# by settgin `error_score` to `np.nan` for instance.

# %%
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import reciprocal

param_distributions = {
    "logisticregression__C": reciprocal(0.001, 10),
    "logisticregression__solver": ["liblinear", "lbfgs"],
    "logisticregression__penalty": ["l2", "l1"],
    "columntransformer__cat-preprocessor__drop": [None, "first"]
}

model_random_search = RandomizedSearchCV(
    model, param_distributions=param_distributions,
    n_iter=200, error_score=np.nan, n_jobs=-1)
model_random_search.fit(df_train, target_train)

# %%
column_results = [f"param_{name}"for name in param_distributions.keys()]
column_results += ["mean_test_score", "std_test_score", "rank_test_score"]

cv_results = pd.DataFrame(model_random_search.cv_results_)
cv_results = cv_results[column_results].sort_values(
    "mean_test_score", ascending=False)
cv_results = cv_results.rename(
    columns={"param_logisticregression__C": "C",
             "param_logisticregression__solver": "solver",
             "param_logisticregression__penalty": "penalty",
             "param_columntransformer__cat-preprocessor__drop": "drop",
             "mean_test_score": "mean test accuracy",
             "rank_test_score": "ranking"})

# %%
cv_results["drop"] = cv_results["drop"].fillna("None")
cv_results = cv_results.dropna("index").drop(columns=["solver"])
encoding = {}
for col in cv_results:
    if cv_results[col].dtype.kind == 'O':
        labels, uniques = pd.factorize(cv_results[col])
        cv_results[col] = labels
        encoding[col] = uniques
encoding
# %%

import plotly.express as px

fig = px.parallel_coordinates(
    cv_results.drop(columns=["ranking", "std_test_score"]),
    color="mean test accuracy",
    dimensions=["C", "penalty", "drop",
                "mean test accuracy"],
    color_continuous_scale=px.colors.diverging.Tealrose,
)
fig.show()


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %%
# TODO add intro with objectives

# ## [markdown]
# Let's first load the data as we did in the previous notebook. TODO add link.

# %%
import pandas as pd

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")
# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

target_name = "class"
target = df[target_name].to_numpy()

data = df.drop(columns=[target_name, "fnlwgt"])

# %% [markdown]
# ## Working with categorical variables
#
# As we have seen in the previous section, a numerical variable is a continuous
# quantity represented by a real or integer number. These variables can be
# naturally handled by machine learning algorithms that are typically composed
# of a sequence of arithmetic instructions such as additions and
# multiplications.
#
# In contrast, categorical variables have discrete values typically represented
# by string labels taken from a finite list of possible choices. For instance,
# the variable `native-country` in our dataset is a categorical variable because
# it encodes the data using a finite list of possible countries (along with the
# `?` symbol when this information is missing):

# %%
data["native-country"].value_counts()

# %% [markdown]
# In the remainder of this section, we will present different strategies to
# encode categorical data into numerical data which can be used by a
# machine-learning algorithm.

# %%
data.dtypes

# %%
categorical_columns = [
    c for c in data.columns if data[c].dtype.kind not in ["i", "f"]]
categorical_columns

# %%
data_categorical = data[categorical_columns]
data_categorical.head()

# %%
print(
    f"The datasets is composed of {data_categorical.shape[1]} features"
)

# %% [markdown]
# ### Encoding ordinal categories
#
# The most intuitive strategy is to encode each category with a different
# number. The `OrdinalEncoder` will transform the data in such manner.

# %%
from sklearn.preprocessing import OrdinalEncoder

encoder = OrdinalEncoder()
data_encoded = encoder.fit_transform(data_categorical)

print(
    f"The dataset encoded contains {data_encoded.shape[1]} features")
data_encoded[:5]

# %% [markdown]
# We can see that the categories have been encoded for each feature (column)
# independently. We can also note that the number of features before and after
# the encoding is the same.
#
# However, one has to be careful when using this encoding strategy. Using this
# integer representation can lead the downstream models to make the assumption
# that the categories are ordered: 0 is smaller than 1 which is smaller than 2,
# etc.
#
# By default, `OrdinalEncoder` uses a lexicographical strategy to map string
# category labels to integers. This strategy is completely arbitrary and often be
# meaningless. For instance suppose the dataset has a categorical variable named
# "size" with categories such as "S", "M", "L", "XL". We would like the integer
# representation to respect the meaning of the sizes by mapping them to increasing
# integers such as 0, 1, 2, 3. However lexicographical strategy used by default
# would map the labels "S", "M", "L", "XL" to 2, 1, 0, 3.
#
# The `OrdinalEncoder` class accepts a "categories" constructor argument to pass
# in the correct ordering explicitly.
#
# If a categorical variable does not carry any meaningful order information then
# this encoding might be misleading to downstream statistical models and you might
# consider using one-hot encoding instead (see below).
#
# Note however that the impact of violating this ordering assumption is really
# dependent on the downstream models (for instance linear models are much more
# sensitive than models built from a ensemble of decision trees).
#
# ### Encoding nominal categories (without assuming any order)
#
# `OneHotEncoder` is an alternative encoder that can prevent the dowstream
# models to make a false assumption about the ordering of categories. For a
# given feature, it will create as many new columns as there are possible
# categories. For a given sample, the value of the column corresponding to the
# category will be set to `1` while all the columns of the other categories will
# be set to `0`.

# %%
print(
    f"The dataset is composed of {data_categorical.shape[1]} features"
)
data_categorical.head()

# %%
from sklearn.preprocessing import OneHotEncoder

encoder = OneHotEncoder(sparse=False)
data_encoded = encoder.fit_transform(data_categorical)
print(
    f"The dataset encoded contains {data_encoded.shape[1]} features")
data_encoded

# %% [markdown]
# Let's wrap this numpy array in a dataframe with informative column names as provided by the encoder object:

# %%
columns_encoded = encoder.get_feature_names(data_categorical.columns)
pd.DataFrame(data_encoded, columns=columns_encoded).head()

# %% [markdown]
# Look at how the "workclass" variable of the first 3 records has been encoded and compare this to the original string representation.
#
# The number of features after the encoding is more than 10 times larger than in the
# original data because some variables such as `occupation` and `native-country`
# have many possible categories.
#
# We can now integrate this encoder inside a machine learning pipeline as in the
# case with numerical data: let's train a linear classifier on
# the encoded data and check the performance of this machine learning pipeline
# using cross-validation.

# %%
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

model = make_pipeline(
    OneHotEncoder(handle_unknown='ignore'),
    LogisticRegression(solver='lbfgs', max_iter=1000))
scores = cross_val_score(model, data_categorical, target)
print(f"The different scores obtained are: \n{scores}")

# %%
print(f"The accuracy is: {scores.mean():.3f} +/- {scores.std():.3f}")

# %% [markdown]
# As you can see, this representation of the categorical variables of the data is slightly more predictive of the revenue than the numerical variables that we used previously.

# %% [markdown]
# ## Exercise 1:
#
# - Try to fit a logistic regression model on categorical data transformed by
#   the OrdinalEncoder instead. What do you observe?
#
# Use the dedicated notebook to do this exercise.

# %% [markdown]
# ## Using numerical and categorical variables together
#
# In the previous sections, we saw that we need to treat data differently
# depending on their nature (i.e. numerical or categorical).
#
# Scikit-learn provides a `ColumnTransformer` class which will send
# specific columns to a specific transformer, making it easy to fit a single
# predictive model on a dataset that combines both kinds of variables together
# (heterogeneously typed tabular data).
#
# We can first define the columns depending on their data type:
# * **binary encoding** will be applied to categorical columns with only too
#   possible values (e.g. sex=male or sex=female in this example). Each binary
#   categorical columns will be mapped to one numerical columns with 0 or 1
#   values.
# * **one-hot encoding** will be applied to categorical columns with more that
#   two possible categories. This encoding will create one additional column for
#   each possible categorical value.
# * **numerical scaling** numerical features which will be standardized.
#
#
#
#
#
#
#

# %%
binary_encoding_columns = ['sex']
one_hot_encoding_columns = [
    'workclass', 'education', 'marital-status', 'occupation',
    'relationship', 'race', 'native-country']
scaling_columns = [
    'age', 'education-num', 'hours-per-week', 'capital-gain',
    'capital-loss']

# %% [markdown]
# We can now create our `ColumnTransfomer` by specifying a list of triplet
# (preprocessor name, transformer, columns). Finally, we can define a pipeline
# to stack this "preprocessor" with our classifier (logistic regression).

# %%
from sklearn.preprocessing import StandardScaler
from sklearn.compose import ColumnTransformer

preprocessor = ColumnTransformer([
    ('binary-encoder', OrdinalEncoder(), binary_encoding_columns),
    ('one-hot-encoder', OneHotEncoder(handle_unknown='ignore'),
     one_hot_encoding_columns),
    ('standard-scaler', StandardScaler(), scaling_columns)])
model = make_pipeline(
    preprocessor, LogisticRegression(solver='lbfgs', max_iter=1000))

# %% [markdown]
# The final model is more complex than the previous models but still follows the
# same API:
# - the `fit` method is called to preprocess the data then train the classifier;
# - the `predict` method can make predictions on new data;
# - the `score` method is used to predict on the test data and compare the
#   predictions to the expected test labels to compute the accuracy.

# %%
from sklearn.model_selection import train_test_split

data_train, data_test, target_train, target_test = train_test_split(
    data, target, random_state=42)
model.fit(data_train, target_train)
model.predict(data_test)[:5]

# %%
target_test[:5]

# %%
data_test.head()

# %%
model.score(data_test, target_test)

# %% [markdown]
# This model can also be cross-validated as usual (instead of using a single
# train-test split):

# %%
scores = cross_val_score(model, data, target, cv=5)
print(f"The different scores obtained are: \n{scores}")

# %%
print(f"The accuracy is: {scores.mean():.3f} +- {scores.std():.3f}")

# %% [markdown]
# The compound model has a higher predictive accuracy than the
# two models that used numerical and categorical variables in
# isolation.

# %% [markdown]
# # Fitting a more powerful model
#
# Linear models are very nice because they are usually very cheap to train,
# small to deploy, fast to predict and give a good baseline.
#
# However, it is often useful to check whether more complex models such as an
# ensemble of decision trees can lead to higher predictive performance.
#
# In the following cell we try a scalable implementation of the Gradient Boosting
# Machine algorithm. For this class of models, we know that contrary to linear
# models, it is useless to scale the numerical features and furthermore it is
# both safe and significantly more computationally efficient use an arbitrary
# integer encoding for the categorical variable even if the ordering is
# arbitrary. Therefore we adapt the preprocessing pipeline as follows:

# %%
from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier

# For each categorical column, extract the list of all possible categories
# in some arbritrary order.
categories = [
    data[column].unique() for column in data[categorical_columns]]

preprocessor = ColumnTransformer([
    ('categorical', OrdinalEncoder(categories=categories),
     categorical_columns)], remainder="passthrough")

model = make_pipeline(preprocessor, HistGradientBoostingClassifier())
model.fit(data_train, target_train)
print(model.score(data_test, target_test))

# %% [markdown]
# We can observe that we get significantly higher accuracies with the Gradient
# Boosting model. This is often what we observe whenever the dataset has a large
# number of samples and limited number of informative features (e.g. less than
# 1000) with a mix of numerical and categorical variables.
#
# This explains why Gradient Boosted Machines are very popular among datascience
# practitioners who work with tabular data.
#
#
#
#
#
#
#

# %% [markdown]
# ## Exercise 2:
#
# - Check that scaling the numerical features does not impact the speed or
#   accuracy of `HistGradientBoostingClassifier`
# - Check that one-hot encoding the categorical variable does not improve the
#   accuracy of `HistGradientBoostingClassifier` but slows down the training.
#
# Use the dedicated notebook to do this exercise.


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
# ---

# %% [markdown]
# In this notebook, we will look at the necessary steps required before any machine learning takes place.
# * load the data
# * look at the variables in the dataset, in particular, differentiate
#   between numerical and categorical variables, which need different
#   preprocessing in most machine learning workflows
# * visualize the distribution of the variables to gain some insights into the dataset

# %%
# Inline plots
# %matplotlib inline

# plotting style
import seaborn as sns

# %% [markdown]
# ## Loading the adult census dataset

# %% [markdown]
# We will use data from the "Current Population adult_census" from 1994 that we
# downloaded from [OpenML](http://openml.org/).

# %%
import pandas as pd

adult_census = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")

# Or use the local copy:
# adult_census = pd.read_csv('../datasets/adult-census.csv')

# %% [markdown]
# We can look at the OpenML webpage to know more about this dataset.

# %%
from IPython.display import IFrame
IFrame('https://www.openml.org/d/1590', width=1200, height=600)

# %% [markdown]
# ## Look at the variables in the dataset
# The data are stored in a pandas dataframe.

# %%
adult_census.head()

# %% [markdown]
# The column named **class** is our target variable (i.e., the variable which
# we want to predict). The two possible classes are `<= 50K` (low-revenue) and
# `> 50K` (high-revenue).

# %%
target_column = 'class'
adult_census[target_column].value_counts()

# %% [markdown]
# Note: classes are slightly imbalanced. Class imbalance happens often in
# practice and may need special techniques for machine learning. For example in
# a medical setting, if we are trying to predict whether patients will develop
# a rare disease, there will be a lot more healthy patients than ill patients in
# the dataset.

# %% [markdown]
# The dataset contains both numerical and categorical data. Numerical values
# can take continuous values for example `age`. Categorical values can have a
# finite number of values, for example `native-country`.

# %%
numerical_columns = [
    'age', 'education-num', 'capital-gain', 'capital-loss',
    'hours-per-week']
categorical_columns = [
    'workclass', 'education', 'marital-status', 'occupation',
    'relationship', 'race', 'sex', 'native-country']
all_columns = numerical_columns + categorical_columns + [
    target_column]

adult_census = adult_census[all_columns]

# %% [markdown]
# Note that for simplicity, we have ignored the "fnlwgt" (final weight) column
# that was crafted by the creators of the dataset when sampling the dataset to
# be representative of the full census database.

# %% [markdown]
# ## Inspect the data
# Before building a machine learning model, it is a good idea to look at the
# data:
# * maybe the task you are trying to achieve can be solved without machine
#   learning
# * you need to check that the data you need for your task is indeed present in
# the dataset
# * inspecting the data is a good way to find peculiarities. These can can
#   arise during data collection (for example, malfunctioning sensor or missing
#   values), or from the way the data is processed afterwards (for example capped
#   values).

# %% [markdown]
# Let's look at the distribution of individual variables, to get some insights
# about the data. We can start by plotting histograms, note that this only
# works for numerical variables:

# %%
_ = adult_census.hist(figsize=(20, 10))

# %% [markdown]
# We can already make a few comments about some of the variables:
# * age: there are not that many points for 'age > 70'. The dataset description
# does indicate that retired people have been filtered out (`hours-per-week > 0`).
# * education-num: peak at 10 and 13, hard to tell what it corresponds to
# without looking much further. We'll do that later in this notebook.
# * hours per week peaks at 40, this was very likely the standard number of
# working hours at the time of the data collection
# * most values of capital-gain and capital-loss are close to zero

# %% [markdown]
# For categorical variables, we can look at the distribution of values:

# %%
adult_census['sex'].value_counts()

# %%
adult_census['education'].value_counts()

# %% [markdown]
# `pandas_profiling` is a nice tool for inspecting the data (both numerical and
# categorical variables).

# %%
import pandas_profiling
adult_census.profile_report()

# %% [markdown]
# As noted above, `education-num` distribution has two clear peaks around 10
# and 13. It would be reasonable to expect that `education-num` is the number of
# years of education. Let's look at the relationship between `education` and
# `education-num`.
# %%
pd.crosstab(index=adult_census['education'],
            columns=adult_census['education-num'])

# %% [markdown]
# This shows that education and education-num gives you the same information.
# For example, `education-num=2` is equivalent to `education='1st-4th'`. In
# practice that means we can remove `education-num` without losing information.
# Note that having redundant (or highly correlated) columns can be a problem
# for machine learning algorithms.

# %% [markdown]
# Another way to inspect the data is to do a pairplot and show how each variable
# differs according to our target, `class`. Plots along the diagonal show the
# distribution of individual variables for each `class`. The plots on the
# off-diagonal can reveal interesting interactions between variables.

# %%
n_samples_to_plot = 5000
columns = ['age', 'education-num', 'hours-per-week']
_ = sns.pairplot(data=adult_census[:n_samples_to_plot], vars=columns,
                 hue=target_column, plot_kws={'alpha': 0.2},
                 height=4, diag_kind='hist')

# %%
_ = sns.pairplot(data=adult_census[:n_samples_to_plot], x_vars='age',
                 y_vars='hours-per-week', hue=target_column,
                 markers=['o',
                          'v'], plot_kws={'alpha': 0.2}, height=12)

# %% [markdown]
#
# By looking at the data you could infer some hand-written rules to predict the
# class:
# * if you are young (less than 25 year-old roughly), you are in the `<= 50K` class.
# * if you are old (more than 70 year-old roughly), you are in the `<= 50K` class.
# * if you work part-time (less than 40 hours roughly) you are in the `<= 50K` class.
#
# These hand-written rules could work reasonably well without the need for any
# machine learning. Note however that it is not very easy to create rules for
# the region `40 < hours-per-week < 60` and `30 < age < 70`. We can hope that
# machine learning can help in this region. Also note that visualization can
# help creating hand-written rules but is limited to 2 dimensions (maybe 3
# dimensions), whereas machine learning models can build models in
# high-dimensional spaces.
#
# Another thing worth mentioning in this plot: if you are young (less than 25
# year-old roughly) or old (more than 70 year-old roughly) you tend to work
# less. This is a non-linear relationship between age and hours
# per week. Some machine learning models can only capture linear interactions so
# this may be a factor when deciding which model to chose.
#
# In a machine-learning setting, we will use an algorithm to automatically
# decide what should be the "rules" in order to make predictions on new data.
# We can check which set of simple rules a decision tree would grasp using the
# same data.

# %%
def plot_tree_decision_function(tree, X, y, ax):
    """Plot the different decision rules found by a `DecisionTreeClassifier`.

    Parameters
    ----------
    tree : DecisionTreeClassifier instance
        The decision tree to inspect.
    X : dataframe of shape (n_samples, n_features)
        The data used to train the `tree` estimator.
    y : ndarray of shape (n_samples,)
        The target used to train the `tree` estimator.
    ax : matplotlib axis
        The matplotlib axis where to plot the different decision rules.
    """
    import numpy as np
    from scipy import ndimage

    h = 0.02
    x_min, x_max = 0, 100
    y_min, y_max = 0, 100
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                         np.arange(y_min, y_max, h))

    Z = tree.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]
    Z = Z.reshape(xx.shape)
    faces = tree.tree_.apply(
        np.c_[xx.ravel(), yy.ravel()].astype(np.float32))
    faces = faces.reshape(xx.shape)
    border = ndimage.laplace(faces) != 0
    ax.scatter(X.iloc[:, 0], X.iloc[:, 1],
               c=np.array(['tab:blue',
                           'tab:red'])[y], s=60, alpha=0.7)
    ax.contourf(xx, yy, Z, alpha=.4, cmap='RdBu_r')
    ax.scatter(xx[border], yy[border], marker='.', s=1)
    ax.set_xlabel(X.columns[0])
    ax.set_ylabel(X.columns[1])
    ax.set_xlim([x_min, x_max])
    ax.set_ylim([y_min, y_max])
    sns.despine(offset=10)

# %%
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams.update(mpl.rcParamsDefault)  # reset the plotting style

# %%
from sklearn.preprocessing import LabelEncoder

# select a subset of data
data_subset = adult_census[:n_samples_to_plot]
X = data_subset[["age", "hours-per-week"]]
y = LabelEncoder().fit_transform(
    data_subset[target_column].to_numpy())

# %% [markdown]
# We will create a decision tree which we will keep really simple on purpose.
# We will only allow a maximum of 2 rules to interpret the results.

# %%
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import plot_tree

max_leaf_nodes = 3
tree = DecisionTreeClassifier(max_leaf_nodes=max_leaf_nodes,
                              random_state=0)
tree.fit(X, y)

# %% [markdown]
# We can now check the set of rules learnt by the tree then visually check
# what theses rules look like.

# %%
# plot the tree structure
fig, ax = plt.subplots()
plot_tree(tree, ax=ax)

# plot the decision function learned by the tree
fig, ax = plt.subplots()
plot_tree_decision_function(tree, X, y, ax=ax)

# %% [markdown]
# Allowing only 3 leaves in the tree, we get similar rules to the ones we
# designed by hand:
# * the persons younger than 28.5 year-old will be considered in the class
#   earning `<= 50K`.
# * the persons older than 28.5 and working more than 40.5 hours-per-week
#   will be considered in the class earning `> 50K`, while the persons working
#   below 40.5 hours-per-week, will be considered in the class
#   earning `<= 50K`.

# %% [markdown]
#
# In this notebook we have:
# * loaded the data from a CSV file using `pandas`
# * looked at the kind of variables in the dataset, and differentiated
#   between categorical and numerical variables
# * inspected the data with `pandas`, `seaborn` and `pandas_profiling`. Data inspection
#   can allow you to decide whether using machine learning is appropriate for
#   your data and to highlight potential peculiarities in your data
#
# Key ideas discussed:
# * if your target variable is imbalanced (e.g., you have more samples from one
#   target category than another), you may need special techniques for machine
#   learning
# * having redundant (or highly correlated) columns can be a problem for
#   machine learning algorithms
# * some machine learning models can only capture linear interaction so be
#   aware of non-linear relationships in your data


text:
# ---
# jupyter:
#   jupytext:
#     formats: python_scripts//py:percent,notebooks//ipynb
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.2'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# #  Exercise 01
# The goal is to write an exhaustive search to find the best parameters
# combination maximizing the model performance

# %%
import pandas as pd

from sklearn.model_selection import train_test_split
from sklearn.preprocessing import OrdinalEncoder
from sklearn.model_selection import RandomizedSearchCV
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
# This line is currently required to import HistGradientBoostingClassifier
from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier

from scipy.stats import expon, uniform
from scipy.stats import randint

df = pd.read_csv(
    "https://www.openml.org/data/get_csv/1595261/adult-census.csv")
# Or use the local copy:
# df = pd.read_csv('../datasets/adult-census.csv')

target_name = "class"
target = df[target_name].to_numpy()
data = df.drop(columns=target_name)

df_train, df_test, target_train, target_test = train_test_split(
    data, target, random_state=42)

from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OrdinalEncoder

categorical_columns = [
    'workclass', 'education', 'marital-status', 'occupation',
    'relationship', 'race', 'native-country', 'sex']

categories = [data[column].unique()
              for column in data[categorical_columns]]

categorical_preprocessor = OrdinalEncoder(categories=categories)

preprocessor = ColumnTransformer(
    [('cat-preprocessor', categorical_preprocessor, categorical_columns)],
    remainder='passthrough', sparse_threshold=0)

from sklearn.experimental import enable_hist_gradient_boosting
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.pipeline import make_pipeline

model = make_pipeline(
    preprocessor, HistGradientBoostingClassifier(random_state=42))

# %% [markdown]
# TODO: write your solution here
#
# Use the previously defined model (called `model`) and using two nested `for`
# loops, make a search of the best combinations of the `learning_rate` and
# `max_leaf_nodes` parameters. In this regard, you will need to train and test
# the model by setting the parameters. The evaluation of the model should be
# performed using `cross_val_score`. We can propose to define the following
# parameters search:
# - `learning_rate` for the values 0.01, 0.1, and 1;
# - `max_leaf_nodes` for the values 5, 25, 45.

# %%
from sklearn.model_selection import cross_val_score

learning_rate = [0.01, 0.1, 1, 10]
max_leaf_nodes = [5, 25, 45]

best_score = 0
best_params = {}
for lr in learning_rate:
    for mln in max_leaf_nodes:
        model.set_params(
            histgradientboostingclassifier__learning_rate=lr,
            histgradientboostingclassifier__max_leaf_nodes=mln
        )
        scores = cross_val_score(model, df_train, target_train, cv=3)
        if scores.mean() > best_score:
            best_score = scores.mean()
            best_params = {'learning-rate': lr, 'max leaf nodes': mln}
print(f"The best accuracy obtained is {best_score:.3f}")
print(f"The best parameters found are:\n {best_params}")


text:
import os
import binascii
import sys
import struct
import math
import hashlib

from Crypto.Cipher import AES
from cryptography.hazmat.primitives.serialization import pkcs12
from cryptography.hazmat.primitives.asymmetric import padding

game_key = b""
content_id = b""
original_iv = b""

file_size = 0
total_blocks = 0

def ReadGameKey(filename): # Read from FAKE.RIF
    global game_key
    global content_id

    riffd = open(filename, "rb")
    riffd.seek(0x50, 0)
    content_id = riffd.read(0x24)
    riffd.seek(0x120, 0)
    game_key = riffd.read(0x10)
    riffd.close()

def ReadPublisherKey(p12, khapp): # used for PSM Developer Assistant (debug PSSE)
    global game_key
    print("[*] Reading PRIVATE KEY from "+p12)
    pkcs12_file = open(p12, 'rb').read()
    passphrase = b"password"
    private_key, certificate, additional_certificates = pkcs12.load_key_and_certificates(pkcs12_file, passphrase) # Parse PKCS12 file

    print("[*] Reading HKAPP "+khapp)
    hkfd = open(khapp, "rb")
    offset_by = 0

    magic = hkfd.read(0x4)
    if magic == b"PAKR": # Keyring
        offset_by = 0x48
    elif magic == b"tkdb": # Protected_kconsole.dat- needs more research
        offset_by = 0x10
    elif magic == b"PSHK": # Raw HKAPP
        offset_by = 0
    else:
        print("[*] Invalid khapp, magic: "+magic.decode("UTF-8"))
        exit()

    hkfd.seek(0x200+offset_by, 0)
    enc = hkfd.read(0x100)
    dec = private_key.decrypt(enc, padding.PKCS1v15()) # RSA Decrypt

    key = dec[0xC0:0xD0]
    iv = dec[0xD0:0xE0]

    hkfd.seek(0x80+offset_by, 0)
    ebuffer = hkfd.read(0x240)
    cipher = AES.new(key, AES.MODE_CBC, iv)
    dec = cipher.decrypt(ebuffer) # AES Decrypt

    # Dont know what these are for.
    something1 = dec[0x00:0x10]
    something2 = dec[0xA0:0xB0]

    game_key = dec[0x60:0x70]

    hkfd.close()

def IvRoll(location):
    iv_next = location
    iv = bytearray(original_iv)
    new_iv = bytearray(iv_next.to_bytes(0x10, "little"))
    for i in range(0,0x10):
        new_iv[i] ^= iv[i]

    return bytes(new_iv)

def GetBlock(fd, block_id):
    global game_key
    global file_size
    global total_blocks
    current_iv = IvRoll(block_id)
    block_loc = block_id * 0x8000
    total_read = 0x8000
    trim_to = total_read

    if block_id == 0: # Skip to filedata
        block_loc = 0x680
        total_read -= 0x680
        trim_to = total_read
    elif block_loc % 0x80000 == 0: # Skip signature block
        block_loc += 0x400
        total_read -= 0x400
        trim_to = total_read

    rd_amt = ((block_loc - 0x680) - (0x400*math.floor(block_loc / 0x80000))) # Total amount of file read so far.

    if block_id >= total_blocks-1: # Is this the last block?
        total_read = file_size - rd_amt
        trim_to = total_read
        total_read += (0x10-(total_read % 0x10))

    if block_id > total_blocks-1:
        return b""
    if total_read <= 0:
        return b""

    fd.seek(block_loc, 0)

    file_data = fd.read(total_read)
    cipher = AES.new(game_key, AES.MODE_CBC, current_iv)
    return cipher.decrypt(file_data)[:trim_to]

def DecryptFile(input_filename):
    global game_key
    global content_id
    global original_iv
    global file_size
    global total_blocks

    psm_dev = False

    iv = b"\x00\x01\x02\x03\x04\x05\x06\x07\x08\x09\x0A\x0B\x0C\x0D\x0E\x0F" # Header IV, used to find the real IV
    key = b"\x4E\x29\x8B\x40\xF5\x31\xF4\x69\xD2\x1F\x75\xB1\x33\xC3\x07\xBE" # Header key for PSM Runtime and used in all retail games.
    psm_dev_key = b"\x00\x11\x22\x33\x44\x55\x66\x77\x88\x99\xAA\xBB\xCC\xDD\xEE\xFF" # Header key for PSM Dev, used for all files encrypted with psm_encryptor.dll
    runtime_game_key = b"\xA8\x69\x3C\x4D\xF0\xAE\xED\xBC\x9A\xBF\xD8\x21\x36\x92\x91\x2D" # key for the PlayStation Mobile Runtime (eg used in Sce.PlayStation.Core.dll)

    total_blocks = math.floor(os.path.getsize(input_filename) / 0x8000)+1

    fd = open(input_filename, "rb")

    # Check magic is "PSSE" or "PSME"
    header = fd.read(4)
    if header != b"PSSE" and header != b"PSME":
        print("[*] "+input_filename.decode("UTF-8")+" is not a valid PSSE File")
        print("[*] Is it already decrypted?")
        exit()

    # Check PSSE Version is 0x1 (not sure what it would be otherwise)
    fd.seek(0x4, 0)
    version = struct.unpack('i', fd.read(4))[0]
    if version != 0x01:
        print("[*] Unsupported PSSE version: "+str(version))
        exit()

    fd.seek(0x8, 0)
    file_size = struct.unpack('Q', fd.read(8))[0]

    fd.seek(0x10, 0)
    psse_type = struct.unpack('i', fd.read(4))[0]
    if psse_type != 0x1:
        print("[*] Unknown PSSE Type: "+str(psse_type))
        exit()

    fd.seek(0x14, 0)
    read_content_id = fd.read(0x24)
    if read_content_id == b"IP9100-NPXS10074_00-0000000000000000": # Sce.PlayStation.Core.dll
        game_key = runtime_game_key
    elif read_content_id == (b"\x00"*0x24): # PSM developer assistant doenst include a Content ID
        psm_dev = True
    elif content_id != read_content_id:
        print("[*] Content ID Mismatch! Expected: "+content_id.decode("UTF-8")+" but got "+read_content_id.decode("UTF-8"))
        exit()

    fd.seek(0x40, 0)
    file_md5 = fd.read(0x10)

    # Supposadly this contains the filename.
    #fd.seek(0x50, 0)
    #enc1 = fd.read(0x20)
    #cipher = AES.new(key, AES.MODE_CBC, iv)
    #dec1 = cipher.decrypt(enc1)

    # Determine file IV
    fd.seek(0x70, 0)
    enc = fd.read(0x10)

    if psm_dev:
        key = psm_dev_key

    cipher = AES.new(key, AES.MODE_CBC, iv)
    original_iv = cipher.decrypt(enc)

    total_file_data = b""

    for i in range(0, total_blocks):
        file_block = GetBlock(fd, i)
        total_file_data += file_block

    fd.close()

    total_file_data = total_file_data[:file_size] # trim to file_size

    md5 = hashlib.md5()
    md5.update(total_file_data)
    got_md5 = md5.digest()

    if file_md5 != got_md5: # Verify decrypted file.
        print("[*] MD5 Mismatch, Expected: "+binascii.hexlify(file_md5).decode("UTF-8")+" got: "+binascii.hexlify(got_md5).decode("UTF-8"))
        print("[*] The game key is most likely wrong.")
        exit();

    return total_file_data

if len(sys.argv) <= 1:
    print("PSSE Decryptor by SilicaAndPina v3")
    print("(for retail games)     <PSM_GAME_FOLDER>")
    print("(for PSM Dev Packages) <PSM_GAME_FOLDER> <PSM_PUBLISHER_KEY> <PSM_HKAPP_OR_KRNG>")
    print("(for anything)         <PSM_GAME_FOLDER> <PSM_GAME_KEY>")
    print("")
    print("For retail games, you require the RO/License/FAKE.RIF file")
    print("And for PSM Dev Packages, you need there original publisher key and keyring.")
    exit()

file = sys.argv[1]
fpath = file.encode("UTF-8")

# Some dumb dirbusting shit.
applications_folder = b"/RO/Application/"
license_file = os.path.normpath(fpath+b"/RO/License/FAKE.rif")
psse_list = os.path.normpath(fpath+applications_folder+b"psse.list")

if not os.path.exists(psse_list):
    applications_folder = b"/Application/"
    license_file = os.path.normpath(fpath+b"/License/FAKE.rif")
    psse_list = os.path.normpath(fpath+applications_folder+b"psse.list")

    if not os.path.exists(psse_list):
        applications_folder = b"/"
        license_file = os.path.normpath(fpath+b"/../License/FAKE.rif")
        psse_list = os.path.normpath(fpath+applications_folder+b"psse.list")

        if not os.path.exists(psse_list):
            print("[*] Cannot find psse.list.")
            exit()

# Direct specify direct game key
if len(sys.argv) == 3:
    game_key = binascii.unhexlify(sys.argv[2])
    f = open(psse_list, "rb")
    f.seek(0x14, os.SEEK_SET)
    content_id = f.read(0x24)
    f.close()

# Read from publisher license
if game_key == b"":
    if len(sys.argv) == 4:
        pkcs12name = sys.argv[2]
        khappname = sys.argv[3]
        print("[*] Reading Publisher Key and App Key")

        if not os.path.exists(pkcs12name):
            print("[*] Publisher Key not found")
            exit()

        if not os.path.exists(khappname):
            print("[*] App key file not found.")
            exit()

        ReadPublisherKey(pkcs12name, khappname)

# Read from fake.rif
if game_key == b"":
    if os.path.exists(license_file):
        ReadGameKey(license_file)

 # Is this a particular file rather than a games folder?
if os.path.isfile(fpath):
    print("Decrypting: "+file)
    data = DecryptFile(file)
    open(fpath, "wb").write(data)
    exit()

# If game key is still unknown after all of this, abandon all hope!
if game_key == b"":
    print("[*] Game key is unknown, i cant decrypt this ! :(")
    exit()

print("[*] psse.list: "+psse_list.decode("UTF-8"))
print("[*] Using game key: "+binascii.hexlify(game_key).decode("UTF-8"))

file_data = DecryptFile(psse_list)
file_list = file_data.replace(b"\r", b"").split(b"\n")

for psse_file in file_list:
    if psse_file == b"":
        continue
    file_path = os.path.normpath(fpath+applications_folder+psse_file)
    print("[*] Decrypting: "+file_path.decode("UTF-8"))
    if os.path.exists(file_path):
        decrypted_data = DecryptFile(file_path)
        open(file_path, "wb").write(decrypted_data)
    else:
        print("[*] Error: File not Found: "+file_path.decode("UTF-8"))

open(psse_list, "wb").write(file_data)
print("[*] Decryption Complete!")

repo_name: LiEnby/psse-decrypt

repo_language: C

mime_type: text/x-c


text:
# psse-decrypt
Decrypt PSSE layer of PSM Games (on PC)
Works on Unity and PSM games, and meets all requirements of: https://github.com/vita-nuova/bounties/issues/73
Usage: ``psse_decrypt.py <psm game folder>``

Discoveries:
 - Every game has its own aes-128-cbc key, (located 0x120 bytes into FAKE.rif), so its probably not possible to decrypt a game without a license
 - PSM Runtime has a hard coded key for its runtime libaries. - A8693C4DF0AEEDBC9ABFD8213692912D

After 9 years, the encryption used by PlayStation Mobile games is finally cracked!!

Requires:
  pycryptodome
  cryptography


text:
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>

typedef struct{
	uint64_t size;        // Size of input in bytes
	uint32_t buffer[4];   // Current accumulation of hash
	uint8_t input[64];    // Input to be used in the next step
	uint8_t digest[16];   // Result of algorithm
}MD5Context;

void md5Init(MD5Context *ctx);
void md5Update(MD5Context *ctx, uint8_t *input, size_t input_len);
void md5Finalize(MD5Context *ctx);
void md5Step(uint32_t *buffer, uint32_t *input);

uint8_t* md5String(char *input);
uint8_t* md5File(FILE *file);

uint32_t F(uint32_t X, uint32_t Y, uint32_t Z);
uint32_t G(uint32_t X, uint32_t Y, uint32_t Z);
uint32_t H(uint32_t X, uint32_t Y, uint32_t Z);
uint32_t I(uint32_t X, uint32_t Y, uint32_t Z);

uint32_t rotateLeft(uint32_t x, uint32_t n);