#TODO add to cointainer
!pip install umap-learn

import numpy as np
from numpy.random import seed

import scipy.stats as stats
from scipy.stats import uniform, truncnorm, randint
import random

import pandas as pd



from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split, GridSearchCV, RandomizedSearchCV
from sklearn.metrics import f1_score, classification_report, confusion_matrix, accuracy_score, roc_auc_score
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error, make_scorer
from sklearn.linear_model import LogisticRegression
import sklearn.cluster as cluster
import umap

from joblib import load
from io import BytesIO
import requests

import matplotlib.pyplot as plt
from pylab import rcParams
import seaborn as sns
plt.style.use('seaborn-poster')
%matplotlib inline

import os 
import pickle
import warnings
from time import time
warnings.filterwarnings('ignore') 

"""
Explanation: Exploratory Data Analysis and Feature Engineering Workshop
Created by:

Vladimir Rybakov, Head of Data Science at Wave Access
Pieter Buteneers (@PieterButeneers), Director of Engineering in ML & AI at sinch.com 

In this workshop you will learn about different methods of data preprocessing and analysis.
Agenda:


Cleaning and transforming data, preparing it for further processing: removing missing values, correcting errors, removing outliers and errors, converting data types, etc. 


Exploratory data analysis: correlation analysis, variable relationship analysis. 


Preparing the feature space: encoding categorical variables, manual feature generation, dimensionality reduction and feature selection.


Training the model.


For this workshop we will be using a dataset, which contains records of approximately 10 000 applications scraped from Google Play Store. This  will allow us to demonstrate different data processing and feature engineering techniques in conditions similar to a real project (such, that data scientists have to deal with in their work all the time).
We will try to answer the question: what factors determine the rating of a particular application and whether this rating can be predicted in advance?
<a id='attention'></a>
!! Attention !!
In order to slightly simplify the process, one assumption was made. We will perform all actions on the whole dataset and will split it into training and test subsets only before actually fitting the model. This approach allows us to avoid writing a large number of functions to convert the test dataset to the format of the training and simplifies the understanding of the performed actions. 
However, it is only acceptable on educational tasks. When working with real data, test and validation datasets should be allocated in advance. All data analysis is performed only on the training subset (except checking data distributions). The validation and test subsets are transformed to the training format (to the same format: only those categories that are in the training set are left, missing values are processed in the same way, etc.)
Libraries
We will use the libraries dividing into the following categories:
<a id='lib_data_manipulation'></a>
data preprocessing and descriptive statistics

pandas - a library for data manipulation: dataset transformation, descriptive statistics, working with time series, simple visualization.
numpy - a library for working with numeric arrays, convenient in combination with pandas, as it uses less resources.
scipy - a library that allows applying mathematical functions, linear algebra, statistical tests to data.
statsmodels - an add-on over SciPy that simplifies the application of some statistical models and makes visualization more illustrative.

<a id='lib_data_visualization'></a>
data visualization

matplotlib - the main library for data visualization, has a rich functionality, although interface is a little bit outdated.
seaborn - an add-on over matplotlib, makes graphics more informative and visually attractive.

<a id='lib_machine_learning'></a>
machine learning

scikit-learn (sklearn) - a library containing a large number of algorithms for solving machine learning problems: data preprocessing, implementation of different models, toy datasets.
umap - implementation of the "Uniform manifest Approximation and Projection" algorithm. Allows you to visualize multidimensional data by lowering its dimension
End of explanation
"""


DATA_F ='./data/googleplaystore_alter.csv'
FEATURE_IMPORTANCE_F= './data/feature_importances_logreg.csv'
BONUS_F = './data/bonus_df_alter.csv'
CLUSTER_F = './data/cluster.joblib?raw=true'

"""
Explanation: <a id='1_data_cleaning'></a>
1. Exploring and  transforming the variables
<a id='surface_examination'></a>
1.1. Brief data overview
First, we should load the data and visually explore at least a few lines, see the types of variables, the number of samples in order to get the general understanding of data and get the feeling of what you are working with.
Let's define the paths to the data. It is better to do this at the beginning so, if in the future we need to change them, you won't have to search through the entire script.
End of explanation
"""


df = pd.read_csv(DATA_F)

"""
Explanation: Construction "./" means we use the current folder of the script. "../" would mean - one level higher relative to the script. In our case we want to stay in the same folder where the script file is located, and go to the "data" folder. 
Pandas uses two data structures to store data:
Series - is an object similar to ordinary one-dimensional array, with the difference that each element has its own index. Moreover, Series has a stored data type.
DataFrame - is a tabular data structure where each column is a Series.
Let's load googleplaystore_alter.csv into variable df with pandas function read_csv().
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# df.head(?)
# df.tail(?)

"""
Explanation: Task 1: Bring out the first (or last) ten rows from dataset using methods head() or tail().
End of explanation
"""


df.iloc[2]

"""
Explanation: Dataset Descriptions:


App - application name;


Category - application category;


Rating - application rating;


Size - application size;


Type - application type (paid or free);


Price - application price;


Content Rating - the age group for which the application is intended;


Genres - subcategories to which the application belongs;


Current Ver - current version of application;


Android Ver - the minimum OS version required for the application works.


DataFrame objects have dual indexing - by columns and by rows. 
To select a specific row you can use loc (returns rows by a given index) or iloc (returns rows by their position in the dataset) methods:
End of explanation
"""


df['App'].head(5)

"""
Explanation: It should be noticed that loc and iloc methods are most likely to return different values. For example, let's select the 3rd line.
We can also select a specific column. With the head function you can print the colum names:
End of explanation
"""


df.at[1, 'App']

"""
Explanation: The at method is used for quick selection of the specific element:
End of explanation
"""


df.loc[1:3, ['App', 'Category']]

"""
Explanation: If you need to select several columns, you could do it by feeding a list. 
For several rows you can use a list:
End of explanation
"""


df.iloc[2:4, [1, 2]]

"""
Explanation: Or you can use a slice:
End of explanation
"""


df[df['Category']=='ART_AND_DESIGN'].head(5)

"""
Explanation: If you need strings that match a condition (for example, only applications of the ART_AND_DESIGN category), the following query format is used (it also returns a DataFrame object):
End of explanation
"""


df[(df['Category']=='ART_AND_DESIGN') & (df.Type == 'Free')].head(5)

"""
Explanation: Here we used the head() function again to only show the first 5 rows).
For complex conditions, each expression is placed in parentheses and the & or | operators are used. For example, we can find the applications from the category ART_AND_DESIGN with the type Free:
End of explanation
"""


df.info()

"""
Explanation: Let's take a look at the general information about dataset using the info() method.
End of explanation
"""


df.Category.value_counts()

"""
Explanation: The dataset consists of 13 columns, only one of which is numeric (float64) and the rest are categorical (object). The dataset has 10841 rows. Note that there are also missing values, mostly in the Current Ver.
Pandas contains many built-in functions for working with data. You can learn more about them here: pandas
<a id='data_studying'></a>
1.2. Exploring the data
The data may contain incorrect or missing values and other artifacts. Moreover, types may not match the content. In order to proceed, we must find and correct all errors. In addition, it is necessary to individually explore every features distribution. 
Based on the research results and some logical considerations, we should decide how to use those features: remove, leave, replace continuous values with classes or categories, reduce the existing number of categories, apply a mathematical function to the values, and so on.
<a id='category'></a>
Category
Let's look at the distribution of applications by category using the value_counts() method:
End of explanation
"""


df[df.Category == '2']

"""
Explanation: Most applications are in FAMILY, GAME and TOOLS  categories.
There is only one application in  "2" category. In addition, this name does not make sense in this context. Let's take a look at this app:
End of explanation
"""


df = df[df.Category != '2']

"""
Explanation: There was a shift of a row on one column to the left. Let's delete this sample.
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# df.Category. ?

"""
Explanation: Task 2: Are there apps with missing category values? You can use methods: isnull and  sum
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# df.Rating. ?

plt.figure(figsize=(7,7))
plt.pie(df.Rating.value_counts(),
        labels=df.Rating.value_counts().index,
        autopct='%1.1f%%',
        startangle=120, 
        explode=[0.02]*3)
plt.axis('equal')
plt.show()

# TODO
# df.Rating.astype(int)
# df.Rating.plot.pie()

df.Rating.describe()

"""
Explanation: <a id='rating'></a>
Rating
Rating - is our target column. It's important to check that all the values are alright.
Task 3 What rating does the majority of apps have?
End of explanation
"""


df.Rating.isnull().sum()

"""
Explanation: Checking for missing values:
End of explanation
"""


df.Size.value_counts()

"""
Explanation: <a id='size'></a>
Size
The information contained in Size is clearly numeric, however the column is categorical because it contains special characters.
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# def check_size(size_value):
#  write your code here  

"""
Explanation: The size of applications is specified in megabytes and kilobytes. But the most common is the value 'Varies with device'.
Because of the characters M (megabytes) and K (kilobytes), pandas will not be able to automatically translate them into numbers and will throw an error. So, it's not possible just change the column type "Size" from object to float using the astype('float') method. 
Let's transfrom all values to the amount of megabytes and then to the number format.
Task 4 Write a function that will receive one of the possible values of the Size column to the input, and will return either the corresponding number of megabytes or np.NaN. Use function string.replace(). 
1Mb = 1024Kb
End of explanation
"""


df.Size = df.Size.apply(check_size)

"""
Explanation: Now let's apply our function to the Size column:
End of explanation
"""


assert df.Size.dtype == float
df.Size.dtype

"""
Explanation: Let's check that the column is in the float format now:
End of explanation
"""


df.Size.isnull().sum()/len(df)

"""
Explanation: All values 'Variations with device' are missing now.
End of explanation
"""


df['unknown_size'] = df.Size.isnull()

"""
Explanation: As you can see there is about 16% of the samples with NaN values (previously 'Variations with device'). That's quite a lot. To deal with this, we can create a new column with binary values that will store information about these NaN samples (perhaps this information will be useful) and then process them in the original column.
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# df.Size = df.Size.fillna(value=?)

"""
Explanation: Now let's handle the missing values. There are several traditional methods to deal with them:

remove
replace with average, mode, median, 0
replace with nearby examples (mostly for time series data)
use the distribution of existing data and use it to generate new values for features. (useful when using linear models)
make this column a target variable and train a separate model to predict the missing values

The easiest option is to delete these rows. But 16% is a lot of samples. We don't want to lose so much useful information. The best way is to use a special model, but it can take a lot of time. Usually it's better to compare different approaches.
For this case we will use the median, but you can try experimenting with other options.
Task 5 Fill in the missing values in the Size column with the median value of the same column. Use the methods fillna и median
End of explanation
"""


assert df.Size.isnull().sum() == 0
df.Size.isnull().sum()

"""
Explanation: Just in case, check that there are no empty values:
End of explanation
"""


df.Type.value_counts(normalize=True)

"""
Explanation: <a id='type'></a>
Type
End of explanation
"""


df.Type.isnull().sum()

"""
Explanation: 93% of apps are free. To say more, you need to look at the relationships with other variables. We will deal with this in the next section. Let's check for missing values:
End of explanation
"""


df[df.Type.isnull()]

"""
Explanation: There is one missing, let's inspect this one.
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")

assert df.Type.isnull().sum() == 0
df.Type.isnull().sum()

"""
Explanation: Task 6 Handle the instance with the missing "Type" in the best way on your opinion
End of explanation
"""


df.Price.value_counts()[:10]

"""
Explanation: <a id='price'></a>
Price
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# df['Price'] = df['Price'].str.replace(?)
# df['Price'] = df['Price']. ?

assert df.Price.dtype == float

"""
Explanation: The variable is numeric, but it contains special characters. That's why it has 'object' Type.
Task 7: Remove the special character "$". Change the column type to 'float' using the method astype().
End of explanation
"""


plt.figure(figsize=[15, 7])
sns.distplot(df.Price[df.Type == 'Paid'])
plt.xlabel("Price, $")
plt.title('The distribution of apps by price' ,size = 24)
plt.show()

"""
Explanation: Let's create a distribution plot of prices among the paid apps.
End of explanation
"""


plt.figure(figsize=[15, 7])
sns.distplot(np.log(df.Price[df.Type == 'Paid']))
plt.xlabel("log(Price), $")
plt.title('The distribution of apps by price',size = 24)
plt.show()

"""
Explanation: Cheap apps make up the most part. Moreover, the rest of apps take too small part of all apps, so it is difficult to say anything by the graphic. 
Let's look at logarithm of price.
End of explanation
"""


df[df.Price > 200]

"""
Explanation: Most apps cost around 1 and 3 dollars. The dispersion of values on the left part of the plot is due to the strong discreteness of small values of the original column. Let's have a look at the strange peak on the right side of the plot.
End of explanation
"""


df.Price.isnull().sum()

"""
Explanation: Strange apps for the rich. Perhaps they are an indicator of status in certain circles or something like that.
Let's check for missing values and move on.
End of explanation
"""


df['Content Rating'].value_counts()

"""
Explanation: <a id='content_rating'></a>
Content Rating
End of explanation
"""


df = df[(df['Content Rating'] != 'Adults only 18+') & (df['Content Rating']!= 'Unrated')] 

df['Content Rating'].value_counts()

"""
Explanation: Remove the "Adults only 18+" and "Unrated classes". They are too few.
End of explanation
"""


df['Content Rating'].isnull().sum()

"""
Explanation: Check for missing.
End of explanation
"""


df.Genres.value_counts()

"""
Explanation: <a id='genres'></a>
Genres
End of explanation
"""


list_of_genres = []
for i in df.Genres.str.split(';').values:
    list_of_genres.extend(i)
print('Amount of subcategory: {}'.format(len(set(list_of_genres))))

"""
Explanation: THere are quite a lot of genres - 119. Moreover, some of them are actually a result of combining two basic genres. Let's find out the amount of unique ones.
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# counts_of_genres = df.Genres.str.findall ?

assert sum(counts_of_genres) == 11333

counts_of_genres.describe()

(counts_of_genres>1).sum()/len(counts_of_genres)

"""
Explanation: Task 8 Calculate the amount of genres for each app and estimate the general distribution for them? Use series.str.findall()
Tip: use the apply function to find the length to keep this a one-liner.
End of explanation
"""


df[counts_of_genres>1].head(10)

"""
Explanation: There isn't a lot of applications with many genres. Let's look at some of them.
End of explanation
"""


df.Genres.isnull().sum()

"""
Explanation: You may notice that the genre column often contains the same as the category column. We will explore this in more detail later. Now let's check for missing values and move on.
End of explanation
"""


df['Current Ver'].value_counts()

"""
Explanation: <a id='current_ver'></a>
Current Version
End of explanation
"""


df = df.drop(columns=['Current Ver'])

"""
Explanation: Because of too many unique values, this isn't a particularly informative column. It seems that each developer uses its own notation, so we will just remove it.
End of explanation
"""


df['Android Ver'].value_counts(normalize=True)

df['Android Ver'].isnull().sum()

"""
Explanation: <a id='android_ver'></a>
Android Version
End of explanation
"""


df = df.dropna(subset=['Android Ver'])

"""
Explanation: There are 2 missing values. Remove it.
End of explanation
"""


for i in range(1, 9):
    df.loc[df['Android Ver'].str.contains('^{}..*'.format(i)), 'Android Ver'] = '{} and up'.format(i)

df['Android Ver'].value_counts(normalize=True).sort_index()

"""
Explanation: Let's cluster them together in bigger classes.
End of explanation
"""


df.loc[df['Android Ver'].str.contains('^[123]..*'), 'Android Ver'] = '1 and up'
df.loc[df['Android Ver'].str.contains('^[5678]..*'), 'Android Ver'] = '5 and up'

df['Android Ver'].value_counts(normalize=True)

"""
Explanation: Groups 1-3 and 5-8 still have a small number of examples. The biggest one is 4 group. We probably should cluster the classes even more and combine 1-3 and 5-8 into two separate groups.
End of explanation
"""


df['Android Ver'].isnull().sum()

"""
Explanation: With such distribution of values, this feature can actually contribute to the model.
End of explanation
"""


df.isnull().sum()

"""
Explanation: <a id='final_check'></a>
Final check
Сheck if there are missing values in the dataset.
End of explanation
"""


df_dup = df[df.duplicated(subset='App')]

df_dup.head()

dup_apps = ['Box', 'Call Blocker', 'Bubble Shooter', 'Word Search']

df_tmp = pd.DataFrame()
for col in dup_apps:
    df_tmp = pd.concat((df_tmp, df[df.App == col]))
df_tmp

"""
Explanation: <a id='removing_duplicates'></a>
1.3. Removing duplicates
Now we should make sure there are no applications that occur in the dataset more than once, or in other words - duplicate each other.
End of explanation
"""


df = df.sort_values('Size')

"""
Explanation: Some applications have duplicates with different sizes, so we’ll sort by size
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# df = df.drop_duplicates(?)

assert df.shape == (9652, 10)

df.shape

"""
Explanation: Task 9 Remove duplicates with drop_duplicates() (keep applications with the largest size)
End of explanation
"""


df.head()

"""
Explanation: After filtering, we lost ~11% from the total volume of our dataset. Obviously this is an unpleasant measure, but it has to be done.
End of explanation
"""


corr_df = df.corr()
cmap = sns.diverging_palette(220, 10, as_cmap=True)
mask = np.zeros_like(corr_df, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True
plt.subplots(figsize=[15,10])
plt.title('Correlation matrix')
sns.heatmap(corr_df, mask=mask, cmap=cmap, linewidths=.5, annot=True)
plt.show()

"""
Explanation: <a id='2_data_relations'></a>
2. Exploring the relationship between variables
<a id='numerical_features'></a>
2.1. Numerical features
One of the simple and effective methods to explore the relationship between numerical features is correlation analysis. 

Correlation  - is a statistical relationship between two or more random variables (or variables that can be considered as random with some acceptable degree of accuracy). Changes in one or more values lead to a systematic changing of another values.

A correlation matrix is a matrix that contains the values of correlations between variables. In our case the correlation values between all pairs of numerical features will be obtained.
<img src="./images/Correlation_examples.png">
The image above give you some insight into what the correlation values mean. As you can see this only works in sort of a straight line. Non-linear correlations, as shown in the bottom row, are not captured by this method.
End of explanation
"""


df_log = df.copy()
df_log['log_price'] = np.log1p(df_log.Price)
df_log.unknown_size = df_log.unknown_size.astype(int)

"""
Explanation: There is no correlation between the numeric variables.

As the next step we will use pairplot method for our numeric variables. In order to be able to tell something we shall consider price on a logarithmic scale. Let's make a separate dataframe for this.
End of explanation
"""


ax = sns.pairplot(df_log, hue='Rating', 
                  vars = ['Size', 'log_price', 'unknown_size'],
                  plot_kws = {'alpha': 0.6, 's': 80, 'edgecolor': 'w'},
                  diag_kind='hist', diag_kws = {'edgecolor': 'w', 'alpha': 0.6, 'bins': 10}, 
                  size = 4)
ax.fig.suptitle('Pairplot with grouping by rating', y=1.02, size=18)
plt.show()

"""
Explanation: Pairplots compares the distributions of variables in pairs and allow you to make the most common assumptions, which can then be checked and clarified later.
The distributions of the variables are shown on the diagonal with a color breakdown by the selected types. Relations between two variables are drawn above and below the diagonal: each point is one application in our case.
By default, pairplots use KDE. But it is often built incorrectly - instead of smooth "bells" there are thin high peaks, ridges and other artefacts. Therefore, the using of histograms is more reliable.
End of explanation
"""


df_app = pd.DataFrame({'App': df.App,
                       'Rating': df.Rating})
df_app['App_len'] = df_app.App.apply(len)

plt.figure(figsize=(15, 7))
for i in sorted(df_app.Rating.unique()):
    sns.kdeplot(df_app.App_len[df_app.Rating==i],
                shade=True,
                legend=False,)
plt.legend(labels=sorted(df_app.Rating.unique()))
plt.xlabel("Number of characters in App")
plt.title("The distribution of the Apps name",size = 18)
plt.show()

"""
Explanation: Task 10: Explore the graphics. What conclusions can you make?
Answer We can see from graphics the following:
+ Applications with a small size and "0" rating are more frequent than others. Applications with other sizes have more instances with rating "2" 
+ No apps with big size and high price
+ Apps with "2" rating mostly have low prices. The peak on the 2nd basket probably occurred due to the replacement of missing "size" samples with the median.
+ Applications with "2" rating stand out among the applications whose size depends on the device 
<a id='app_len_&_rating'></a>
2.2. Application length and rating
Let's see how the length of the application name affects the rating. To do this, create a separate dataframe with the App and Rating columns. As a result we will receive a column with the length of the application name.
End of explanation
"""


df.groupby('Type')['Price'].mean()

"""
Explanation: Apps with "0" rating are more likely to have a short title. At the same time, more popular applications are in the range of 20 to 50 characters.
<a id='type_&_price'></a>
2.3. Price and Type
It is logical to assume that if "Type" is "Free", then "Price" should be equal to zero.
Let's group "Type" and look at the average price in each group ("Free" and "Paid").
End of explanation
"""


df.head(5)

"""
Explanation: The average price of a "Free" app is zero dollars, a paid one is 14 dollars. Looks fine.
<a id='category_&_genres'></a>
2.4. Category and Genres
Genres is a subsection of Category.
End of explanation
"""


df_cat_genrs = df[['Category', 'Genres', 'Rating']]
df_cat_genrs['Count_of_genres'] = df_cat_genrs.Genres.str.findall(';').apply(len)+1
df_cat_genrs.head()

"""
Explanation: Let's see which categories use the most subcategories.
Next we will aggregate categories and genres in a separate table and make a new column with the number of genres.
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# group_df = df_cat_genrs.groupby( ? )[ ? ].describe().sort_values('mean', ascending=False)
group_df.head(15)

"""
Explanation: Task 11 Explore the genre number statistics for each category. Use groupby и describe.
End of explanation
"""


df_cat_genrs.Category = df_cat_genrs.Category.str.lower()
df_cat_genrs.Genres = df_cat_genrs.Genres.str.lower()
df_cat_genrs.Genres = df_cat_genrs.Genres.str.replace('&', '')
df_cat_genrs.Genres = df_cat_genrs.Genres.str.replace(';', ' ')
df_cat_genrs.Category = df_cat_genrs.Category.str.replace('_', ' ')
df_cat_genrs.Category = df_cat_genrs.Category.str.replace('and', '')

df_cat_genrs.head(5)

"""
Explanation: The maximum number of subcategories in one category is 2, the minimum is 1. Two subcategories are most often found in the categories PARENTING, FAMILY, EDUCATION. Only 11 categories have multiple genres.
Categories and genres often overlap. Let's study this more in detail.
For starters we should bring categories and genres to a common format: transform them into lower case (str.lower()) and remove special characters (replace the special characters with spaces or nothing using the str.replace()).
Note: pandas has many methods for working with strings. You can read more about them in the documentation
Working with Text Data.
End of explanation
"""


(df_cat_genrs['Category'] == df_cat_genrs['Genres']).sum()/len(df_cat_genrs)

"""
Explanation: Let's calculate how many apps have the same genres and categories.
End of explanation
"""


df_cat_genrs['is_cat_equal_genre'] = df_cat_genrs['Category'] == df_cat_genrs['Genres']
df_cat_genrs.head()

"""
Explanation: The columns are completely the same for almost 70%. That is, 70% of the values don't contain new information. Most likely, there are more matches considering possible errors related to usage of regular expressions on raw data.
Let's add a column to indicate whether the genre is equal to the category or not.
End of explanation
"""


plt.figure(figsize=(10, 5))

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# sns.countplot(x= ?,  hue='Rating', data=df_cat_genrs, hue_order=np.sort(df['Rating'].unique()))

plt.xlabel("Coincidence of Genre and Category")
plt.title('The distribution of apps by ratings',size = 18)
plt.show()

"""
Explanation: Task 12 Explore whether there is a difference between distributions of apps with the same and with different categories and genres (use sns.countplot).
End of explanation
"""


cr_tab = pd.crosstab(df_cat_genrs.Rating, df_cat_genrs.is_cat_equal_genre)
cr_tab

cr_tab / cr_tab.min()

"""
Explanation: The rating distributions for True and False are slightly different. You can see that applications with a rating 0 stand out a little bit. To numerically estimate this distribution, we will use the contingency table.
End of explanation
"""


df = df.drop(columns=['Genres'])
df['is_cat_equal_genre'] = df_cat_genrs['is_cat_equal_genre'] 

"""
Explanation: The difference in distribution for apps with same category and genre is negligible. On the other side, applications with different genre and category are 1.5 times more likely to have "2" rating.
Such column could be useful in predicting the rating of applications. At the same time there is little sense to leave two columns with almost the same information. That's why we remove the column genre and add column is_cat_equal_genre in the main dataset.
End of explanation
"""


# for task version - hide next 3 cells

plt.figure(figsize=(10, 5))

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# sns.countplot(x= ?, hue= ?, data=df, hue_order=np.sort(df['Rating'].unique()))

plt.xlabel("Сontent Rating")
plt.title('The distribution of apps by ratings',size = 18)
plt.show()

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# df_ct_rc = pd.crosstab(?, ?)

df_ct_rc /= df_ct_rc.sum()
df_ct_rc

plt.figure(figsize=(10,10))
for i, con_rating in enumerate(df['Content Rating'].unique()):
    plt.subplot(2, 2, i+1)
    vc = df.Rating[df['Content Rating']== con_rating].value_counts().sort_index()
    patches = plt.pie(vc,  autopct='%1.1f%%', 
           startangle=120, explode=[0.02]*3)
    plt.title(con_rating)
    plt.axis('equal')
plt.legend(labels=vc.index, loc=(1.02,1.75))

   
plt.show()

"""
Explanation: <a id='rating_&_content_rating'></a>
2.5. Rating and Content Rating
Task 13 Explore the relationship between rating and content rating columns.
End of explanation
"""


pd_ct_comp = pd.crosstab(df.Rating, [df['Content Rating'], df.is_cat_equal_genre])
pd_ct_comp /= pd_ct_comp.sum()
pd_ct_comp = pd_ct_comp.style.background_gradient(cmap='summer_r')
pd_ct_comp

"""
Explanation: Let's look at the table of both at the rating contingency with Content Rating column and the is_not_equal_genre column.
End of explanation
"""


g = sns.catplot(x="Category",y="log_price", data=df_log[df_log.Type=='Paid'])
g.fig.set_figheight(7)
g.fig.set_figwidth(20)
plt.title('The distribution of logarithmic price in different categories', size = 20)
plt.xticks(rotation=90)
plt.show()

"""
Explanation: First thing you could notice are applications with the age rating "Everyone" and the same categories and genres. For those application the most common rating is "0". Apps with Everyone 10+ and is_not_equal_genre False  mostly have rating "2".
The most interesting thing is that the distribution of the rating for Teen and Mature 17+ applications  are fairly the same in general, they have a different distribution for applications with the same genre and category (is_cat_equal_genre True).
<a id='category_&_price'></a>
2.6. Category and Price
Another convenient way to explore the data distribution is a catplot. In its simplest, it allows you to display each individual example in an individual column. This way you can get a sense of the data density. 
We will use the previously created logarithmic price column and consider only paid applications and categories with at least 10 such applications.
End of explanation
"""


g = sns.catplot(x="Category",y="log_price", data=df_log[df_log.Type=='Paid'], kind='bar')
g.fig.set_figheight(7)
g.fig.set_figwidth(20)
plt.title('The distribution of logarithmic price in different categories', size = 20)
plt.xticks(rotation=90)
plt.show()

"""
Explanation: There are not a lot of paid applications, but numerous categories, so the graphics are so sparse. You can see that some categories have outlier applications: for instance, Communication category applications are mostly grouped together but one stands out from the others. Finance and Lifestyle are also quite sparse. Note, that the Family group contains a group of detached expensive applications.
Now will take a look at the same distributions, but using confidence intervals.

Confidence Intervals - is a type of interval estimate used in statistics that are calculated for a given level of significance. They allow us to make a statement that the true value of an unknown statistical parameter of the general population is in the obtained range of values with a probability that is given by the selected level of statistical significance.
What is the practical meaning of the confidence interval?


A wide confidence interval indicates that the sample average doesn't accurately approximate the general average. This is usually related to an insufficient sample size, or to its heterogeneity, i.e. the large variance. Both give a large average error and, correspondingly, a wide confidence interval. This is the reason for returning to the planning stage of the study. 


The upper and lower limits of the confidence interval give an estimation of whether the results will be statistically significant or not.
End of explanation
"""


df_log.log_price[(df_log.Category.isin(['LIFESTYLE', 'FINANCE'])) & (df_log.Type == 'Paid')].describe()

df_log[(df_log.log_price > 5.) & (df_log.Category.isin(['LIFESTYLE', 'FINANCE']))]

"""
Explanation: In this graph, the column height shows the average value, and the bar shows the confidence interval. You can notice that the Finance and Lifestyle categories have a very wide confidence interval. This means that these categories are highly sparse and have some isolated groups in different parts of the distributions. It is not valid to use the average value for such categories. There is no sense to consider other categories with a wide confidence interval - there are too few examples to interpret the stats meaningfully.
The Family category has a narrow confidence interval: although there is a separate group of expensive applications, their number is small and the distribution of prices in this category is uniform in general.
Let's look at the groups of expensive applications in the categories Finance and Lifestyle
End of explanation
"""


g = sns.catplot(x="Category",y="log_price", data=df_log[(df_log.Type=='Paid') & (df_log.log_price < 5)], kind='bar')
g.fig.set_figheight(7)
g.fig.set_figwidth(20)
plt.title('The distribution of logarithmic price in different categories', size = 20)
plt.xticks(rotation=90)
plt.show()

"""
Explanation: The shift was due to a few expensive "freaky" apps we already saw earlier. Let's remove them and rearrange the graphic.
End of explanation
"""


df = df[~((df_log.Category.isin(['LIFESTYLE', 'FINANCE'])) & (df_log.log_price > 5))]
df_log = df_log[~((df_log.Category.isin(['LIFESTYLE', 'FINANCE'])) & (df_log.log_price > 5))]

"""
Explanation: The Lifestyle and Finance categories became more realistic and stable. So, we are removing extra expensive samples.
End of explanation
"""


df_log['log_size'] = df_log.Size.apply(np.log1p)

"""
Explanation: <a id='Категории_и_размер'></a>
2.7. Category and Size
Let's add the logarithm of size to df_log.
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# g = sns.catplot(x= ?, y= ?, data=df_log)

g.fig.set_figheight(7)
g.fig.set_figwidth(20)
plt.title('The distribution of logarithmic size in different categories', size = 20)
plt.xticks(rotation=90)
plt.show()

"""
Explanation: Task 14 Explore the distribution of app size in each category. Use the df_log table and the log_size column
End of explanation
"""


results = pd.DataFrame(columns=['method', 'model', 'val score', 'test score', 'learning time', 'predict time'])

"""
Explanation: <a id='3_feature_space'></a>
3. Feature Engineering
Let's create a table where we will store the results.
End of explanation
"""


class Timer(object):
    def __init__(self):
        self.elapsed_time = 0
    def __enter__(self):
        self.start = time()
    def __exit__(self, type, value, traceback):
        self.end = time()
        self.elapsed_time = int((self.end - self.start)*1000)

timer = Timer()

"""
Explanation: To measure processing time we will use a special context manager.
End of explanation
"""


Y = df.Rating.astype(int)
df = df.drop(columns=['Rating'])

"""
Explanation: <a id='base_model'></a>
3.1. Base model
Our task is to determine the rating of the application. Let's take this column into a separate variable, and remove it from the main dataset.
End of explanation
"""


y_train, y_test = train_test_split(Y,  test_size=0.3, random_state=42)

"""
Explanation: The so-called majority classifier can be the simplest type of a base model. The point is: we always predict the class that is most common in our training dataset.
Let's split the data into training and test sets. For appropriate comparision of the different models, this partition must always be the same, that's why we use fixed random_state. For now, we don't need the features, only the target variable.
We remind you that the train/test split of dataset samples should actually be performed before data processing!!!
End of explanation
"""


classes, counts = np.unique(y_train, return_counts=True)
major_class = classes[np.argmax(counts)]
major_class

"""
Explanation: Now let's determine which class is most common in the training data.
End of explanation
"""


base_predict = np.full(y_test.shape, major_class)

"""
Explanation: In the training set, samples with second class are more common. Then our test prediction will be an array with the same size as the size of out test set and it will be completely filled with 2.
End of explanation
"""


def calculate_metrics(y_predict, y_test):
    print('acсuracy: {:.4f}'.format(accuracy_score(y_predict, y_test)))
    print('F1 score: {:.4f}'.format(f1_score(y_predict, y_test, average='macro')))
    print(classification_report(y_predict, y_test))
    print(confusion_matrix(y_predict, y_test))

"""
Explanation: To analyze the results, we will use a function that displays several different metrics.
End of explanation
"""


calculate_metrics(base_predict, y_test)

"""
Explanation: Now let's calculate the accuracy of the majority classifier.
End of explanation
"""


results = results.append({'method':'Baseline',
                          'model': 'Majority',
                          'val score': None,
                          'test score': f1_score(base_predict, y_test, average='macro'), 
                          'learning time' : None,
                          'predict time': None},
                          ignore_index=True)

"""
Explanation: These will be our reference values. Let's write them in the table.
End of explanation
"""


df.info()

numeric_df = df.select_dtypes(include=['int64', 'float64', 'bool'])
numeric_df.head()

"""
Explanation: <a id='origin_features'></a>
3.2. Basic features
Now we are going to use the original feature space to predict the target variable and measure the accuracy of the prediction. We will use only numeric and boolean features: they don't require the preparation, unlike text.
End of explanation
"""


X_train, X_test, y_train, y_test = train_test_split(numeric_df, Y, 
                                                    test_size=0.3, 
                                                    random_state=42)

"""
Explanation: Now we'll get our test and training sets. (notice that we use the same random_state and therefore the data is split the same way as before)
End of explanation
"""


pipe = Pipeline([('scale', StandardScaler()),
                 ('clf', LogisticRegression(random_state=42))])

params = {
    'clf__C': [0.01, 0.05, 0.1, 0.5, 0.9, 0.99],
    'clf__penalty': ['l1', 'l2', 'elasticnet']
} 

np.random.seed(123)

clf = GridSearchCV(pipe, 
                   cv=3,  
                   param_grid=params, 
                   scoring='f1_macro',
                   verbose=1,
                   n_jobs=6)

with timer:
    clf.fit(X_train, y_train)

learning_time = timer.elapsed_time

best_clf = clf.best_estimator_.steps[1][1]
best_clf

clf.best_score_

with timer:
    predict = clf.predict(X_test)
    
predict_time = timer.elapsed_time

calculate_metrics(predict, y_test)

"""
Explanation: We will train a Logistic Regression model - one of the simplest ones among linear classifiers. To perform hyperparameter optimization we'll use cross-validation by applying Pipeline and GridSearchCV functions. This specific classifier was chosen in order to decrease the training time within the workshop. Results for more complex models will be provided in the bonus file.

Data scaling — the process of adjusting the data to a single scale. The main methods are normalization - adjusting all features to a value in the range from 0 to 1, and standardization - data preprocessing, after which each feature has an average of 0 and a variance of 1

Scaling the data can significantly affect the process of training a model. For example, it increases the speed and the stability of gradient descent algorithm. Moreover, methods that rely on working in a multidimensional space (such as KNN) will not work correctly if the features have a different scale. These effects don't always appear, but scaling the data is unlikely to make the results worse, so you should always perform it before training the model.
Now let's create a pipeline and set a grid of the hyperparameters for the classifier.
Note: if the Logistic Regression algorithm does not converge because the of the parameters selected it will throw warnings.
End of explanation
"""


results = results.append({'method':'Numeric',
                          'model': 'LR',
                          'val score': clf.best_score_,
                          'test score': f1_score(clf.predict(X_test), y_test, average='macro'), 
                          'learning time': learning_time,
                          'predict time': predict_time},
                          ignore_index=True)

"""
Explanation: Again, we shall save the results in the table.
End of explanation
"""


df_dummies = df.copy()

"""
Explanation: As you can see, the F1-score of the model is slightly higher than the score of the majority classifier.
As the next step we should add categorical features. First we will encode them features into numerical ones. This can be done in multiple ways.
For example, you could take one column, index all unique category names with numbers and just replace them in the column. This type of encoding is called Label Encoding. This method is usually bad because most of the time the values of numbers don't reflect the category relation. For example, encoding the application with type GAME as 1 and with MEDICAL as 2 turns out that MEDICAL > GAME, while such a comparison is incorrect. Therefore, many methods, especially linear ones, will not work well with this type of encoding. At the same time, methods based on decision trees will work fine.
Another popular type of encoding is called Dummy Encoding or One-hot Encoding. Each unique category name has its own column. The values in this column can only be 0 and 1. If the sample instance has some category, then the corresponding dummy-column will contain a "1" value, and all others will be filled with zeros. This type of encoding has the disadvantage that if there are a lot of unique categories, the size of the feature matrix becomes too large and sparse. Therefore additional preprocessing of categorical data may be necessary.
There are other, less popular and well known techniques, but here we will stick to One-hot encoding.
First of all, we need to copy all the data to the new dataset.
End of explanation
"""


df_dummies = df_dummies.drop(columns=['App'])

"""
Explanation: Next, it is necessary to delete the application names - they are unique for each sample. If we encode them with dummy encoding, it will lead to adding N, where N is the length of the whole dataset, and all these columns will have only one value 1 and all the others 0. Such features obviously do not work.
End of explanation
"""


con_rat_dict =  {'Everyone':0, 'Everyone 10+':1, 'Teen':2, 'Mature 17+':3}
df_dummies = df_dummies.replace({"Content Rating": con_rat_dict})

"""
Explanation: As was said before, sometimes it makes sense to encode columns with LabelEncoding. In our case the Content Rating column is suitable for this approach. We should encode it in a way, so with the growth of restrictions, the corresponding number will also grow.
End of explanation
"""


df_dummies = pd.get_dummies(df_dummies)

df_dummies.shape

"""
Explanation: Now let's encode the remaining categories with the dummy method.
End of explanation
"""


X_train, X_test, y_train, y_test = train_test_split(df_dummies, Y,
                                                    test_size=0.3, 
                                                    random_state=42)

"""
Explanation: Now we have 44 features instead of 4.
Let's prepare a training and test set.
End of explanation
"""


np.random.seed(123)

with timer:
##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
    # clf = GridSearchCV(pipe, cv=3,  
                       # param_grid=params, 
                       # scoring='f1_macro', verbose=1,
                       # n_jobs=6).fit(?, ?)
    
learning_time = timer.elapsed_time

print('Best score:', clf.best_score_)

with timer:
    predict = clf.predict(X_test)
    
predict_time = timer.elapsed_time

calculate_metrics(predict, y_test)

results = results.append({'method': 'Categories',
                          'model': 'LR',
                          'val score': clf.best_score_,
                          'test score': f1_score(clf.predict(X_test), y_test, average='macro'), 
                          'learning time': learning_time,
                          'predict time': predict_time},
                          ignore_index=True)

"""
Explanation: Task 15 Train the model, make a prediction on the test data and output the statistics. Don't forget to measure your training time.
End of explanation
"""


df_new = df.copy()

"""
Explanation: As you can see, adding categorical features have improved the accuracy of the model.
<a id='feature_generation'></a>
3.3. Feature generation
Now we will expand our feature space by adding new features - manually crafted with love and care! 
First we need a separate copy of the original dataset.
End of explanation
"""


df_new['price_for_mb'] = df_new.Size/df_new.Price
df_new['price_for_mb'] = df_new['price_for_mb'].replace([np.inf, -np.inf], 0)

"""
Explanation: Let's add a price per megabyte column.
End of explanation
"""


df_new['log_price'] = df_new.Price.apply(np.log1p)
df_new['log_size'] = df_new.Size.apply(np.log1p)

"""
Explanation: Also let's add logarithms for price and size.
End of explanation
"""


p = PolynomialFeatures(degree=2).fit(df_new[['log_price', 'log_size']])
poly_df = pd.DataFrame(p.transform(df_new[['log_price', 'log_size']]), 
                        columns=p.get_feature_names(['log_price', 'log_size']))
poly_df.head(10)

"""
Explanation: Usually it is useful to indroduce polynomial features for the numerical ones. Let's make them from the previously calculated logarithms of price and size. To do this, we can use the scikit-learn function Polynomial Features.
End of explanation
"""


poly_df = poly_df.drop(columns=['log_price', 'log_size', '1'])
poly_df.index = df_new.index
df_new = pd.concat([df_new, poly_df.reindex(df_new.index)], axis=1)

df_new.head(5)

"""
Explanation: You can see that there are new columns low_price^2, log_price, log_size and log_size^2. Now you need to attach them to the main dataset, and discard unnecessary ones. To avoid errors during concatenation, we will replace the poly_df index with the df_new index.
End of explanation
"""


df_new['len_of_app_title'] = df_new.App.apply(len)
df_new['count_of_app_title'] = df_new.App.str.split(' ').apply(len)

"""
Explanation: Add the number of characters and the number of words in the title as features.
End of explanation
"""


# Let's change to lower case, remove special characters, leave only Latin letters and numbers.
df_new['cleantext'] = df_new.App.str.lower()
df_new['cleantext'] = df_new.cleantext.str.replace('[-_]', ' ')
df_new['cleantext'] = df_new.cleantext.str.replace('[^0-9A-Za-z ]+', '')

"""
Explanation: The names of applications contain a lot of garbage. It is better to clean it.
End of explanation
"""


df_new[['App', 'cleantext']].head(10)

"""
Explanation: First we need to check what the text looks like after cleaning.
End of explanation
"""


df_new['len_of_cleantext_title'] = df_new.cleantext.apply(len)
df_new['count_of_cleantext_title'] = df_new.cleantext.str.split(' ').apply(len)

"""
Explanation: Now, after filtering, we can add the number of characters and the number of words as features.
End of explanation
"""


df_new['diff_len_title'] = df_new['len_of_app_title'] - df_new['len_of_cleantext_title']
df_new['diff_count_title'] = df_new['count_of_app_title'] - df_new['count_of_cleantext_title']

df_new[['diff_len_title', 'diff_count_title']].describe()

"""
Explanation: Next step is to add the difference between the original number of words and symbols and the normalized ones.
End of explanation
"""


df_new[['App', 'cleantext']][df_new['diff_count_title'] == -5]

"""
Explanation: For some reason there are normalized headers that have more words than original ones. We need to check it.
End of explanation
"""


len_vocab = len(set(' '.join(df_new.cleantext.tolist()).split(' ')))
max_count_of_cleantext_title = round(df_new.count_of_cleantext_title.max(), 0)
print('Number of unique words: {}'.format(len_vocab))
print('Maximum number of words in a title: {}'.format(max_count_of_cleantext_title))

df_new = df_new.drop(columns=['App', 'cleantext'])

"""
Explanation: Everything seems fine.
End of explanation
"""


df_new = df_new.replace({"Content Rating": con_rat_dict})

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# df_new = pd. ?

df_new.shape

"""
Explanation: Let's encode categorical features into numerical vectors. 
Task 16: Encode categorical features to one-hot format. Don't forget about the special case with Content Rating.
End of explanation
"""


scaled_data = StandardScaler().fit_transform(df_new)

%%time
np.random.seed(1)
X_umap = umap.UMAP(n_components=2, random_state=2).fit_transform(scaled_data)

plt.figure(figsize=[15, 9])
plt.title('UMAP')

for i in Y.unique():
    indx = np.where(Y.values == i)
    plt.scatter(X_umap[indx, 0], X_umap[indx, 1], marker='.', alpha=0.5, label=i)
    plt.legend()
plt.show()

"""
Explanation: Now it's time to visualize our data. We will use UMAP for this.
It is advisable to standardize the data before using this algorithm.

UMAP (Uniform Manifold Approximation and Projection) - this is a relatively new dimensionality reduction algorithm. Previously the most popular algorithm for these purposes was t-SNE. However UMAP is better in almost every way: UMAP has no restrictions on the dimension of the original feature space, it is much faster and more computationally efficient than t-SNE, and also better in the task of transferring the global data structure to a new, reduced space (as the authors of the algorithm claim). For these reasons t-SNE is hardly used these days.
End of explanation
"""


cls = cluster.KMeans(n_clusters=30, random_state=42)
kmeans_labels = cls.fit_predict(scaled_data)

# For Google Colab below

# Unfortunately, Google collab doesn't allow to reproduce KMeans model training in spite of the random_state.
# Without it we can't reproduce our clusters and distributions inside them. That's why we're loading pretrained cluster
# model. If you want train your own cluster model you can use the commented code

# m_file = BytesIO(requests.get(CLUSTER_F).content)
# cls = load(m_file)
# kmeans_labels = cls.predict(scaled_data)

"""
Explanation: You can see that there is a lot of clusters and they contain different proportions of applications with different ratings. Let's apply a clustering algorithm to separate these groups. For now we will assume that we have 30 clusters and use the k_means algorithm.
End of explanation
"""


plt.figure(figsize=[15, 9])
ax = plt.subplot()

for i in np.unique(kmeans_labels):
    indx = np.where(kmeans_labels == i)
    plt.scatter(X_umap[indx, 0], X_umap[indx, 1], marker='.', alpha=0.5, label=i)
ax.legend(loc='upper center', bbox_to_anchor=(1.05, 1.0),  ncol=1, fancybox=True, shadow=True, fontsize=10)
plt.show()

"""
Explanation: Now we visualize the results of our clustering.
End of explanation
"""


df_clstr = pd.DataFrame({'Cluster': kmeans_labels, 'Rating':Y})
g = sns.catplot("Rating", col="Cluster", col_wrap=5, data=df_clstr, kind="count", height=2.5, aspect=.8, size=4)

"""
Explanation: Looks promissing. Let's explore whether there is a difference in the distributions of the target variable in each cluster.
End of explanation
"""


df_clstr = df_clstr.drop(columns=['Rating'])
df_clstr['Cluster'] = df_clstr['Cluster'].astype('category')
df_clstr = pd.get_dummies(df_clstr)
df_clstr.index = df_new.index

df_new=pd.concat([df_new, df_clstr], axis=1)
df_new.shape

"""
Explanation: Obviously, the distribution of the target variable differs in many way throughout the clusters. We probably should add the cluster indexes to the main dataset. Note that the cluster index is a category. To perform one-hot encoding we will remove the target variable and then we can add new columns to the main dataset.
End of explanation
"""


X_train, X_test, y_train, y_test = train_test_split(df_new, Y,
                                                    test_size=0.3, 
                                                    random_state=42)

np.random.seed(123)

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# with ?:
#    clf = GridSearchCV( ? ).fit( ? )

learning_time = timer.elapsed_time

print('Best score:', clf.best_score_)

print('Best params:', clf.best_params_)

with timer:
    predict = clf.predict(X_test)
    
predict_time = timer.elapsed_time

calculate_metrics(predict, y_test)

results = results.append({'method': 'Manual features',
                          'model': 'LR',
                          'val score': clf.best_score_,
                          'test score': f1_score(clf.predict(X_test), y_test, average='macro'), 
                          'learning time':learning_time,
                          'predict time': predict_time},
                          ignore_index=True)

"""
Explanation: It also would be possible to use aggregation features: for example, average, median, minimum, etc. for each category, age rating or other categorical features. You can apply simple mathematical operations to pairs of numerical features (for example, multiply the size of the application by the price): sometimes they don't make explicit sense, but they might be useful. There are also cases of logical addition of features (there is a category Medicine and an age rating of 18+, for example). To obtain such features, it is convenient to use the featuretools library. In general, the number of possible new features is just limited only by your imagination. Unfortunately today we don't have enough time go all crazy with manually creating these feature, so we will limit ourselves to the ones that we have already discussed.
Now let's try see if our work paid off
Task 17: Split the dataset and target variable into training and test set, train the model, make a prediction on the test data and output the statistics. Don't forget to measure your training time.
End of explanation
"""


model = PCA()
pca_data = model.fit_transform(StandardScaler().fit_transform(df_new))

plt.figure(figsize=[13, 5])
plt.title('Principal Component Analysis (PCA)')
plt.plot(range(len(model.explained_variance_ratio_)), model.explained_variance_ratio_, '--o')
plt.ylabel('Explained variance ratio')
plt.xlabel('Components')
plt.tight_layout()
plt.show()

"""
Explanation: Adding the manually created features increased the effectiveness of the model.
<a id='dimensional_reduction'></a>
4. Feature Selection and Dimensionality Reduction
Usualy it is useful to know which feature are actually doing the work, so we can drop all others as they can actually cause over-fitting and thus hurt the final score. In our case there are only 123 features, but it can be tens and hundreds of thousands. The goal is to reduce the dimensionality of the feature space and not lose useful information.
There are various ways to reduce the dimension of the data:

Statistical methods of dimension reduction
1. Principal Component Analysis
2. Linear Discriminant Analysis
3. Factor analysis
*Machine Learning based algorithms
1. T-SNE 
2. UMAP and others.
Methods based on feature selection
1. The significance of the features (determined by the training of some models)
2. Filtering methods (rank features based on statistic functions - e.g. Chi square test, ANOVA)
3. Lasso and Ridge regression
4. Greedy algorithm - step by step adding features (Forward Feature Selection) or Backward Feature Elimination
Autoencoders should be noted as a method of reducing the dimension of data but will not be covered here. 

We will look at a couple of these methods: Principal Component Analysis and Backward Feature Elimination.
<a id='pca'></a>
4.1. Principal Component Analysis
Principal Component Analysis (PCA) - one of the popular ways to reduce data dimensionality with minimum loss of information. The challenge of principal component analysis is to find subspaces of smaller dimension, in orthogonal projection on which the data variation (that is, the standard deviation from the mean value) is maximal.
In simple terms, the method finds the n-th number of new axes, projections on which retain as much variance of the original data as possible. The projections themselves become new features.
It should be used carefully, because the feature space is radically changing. The information necessary for the algorithms performance can be lost even with a high retention of information content. Also it is very prone to overfitting.
Since the method is related to variance of features, it is necessary to scale data before using PCA.
End of explanation
"""


##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# model = PCA(n_components=?)
# pca_data = model.fit_transform(?)

assert pca_data.shape[1] == 48
pca_data.shape

"""
Explanation: Apparently, most of the information can be saved with 60-70 new projections. Let's see how many components do we need to keep 99% of the original dataset's information.
Task 18 Apply the PCA method to the df_new dataset again, but set the n_components parameter to 0.99. If this parameter is an integer, then in the final matrix will be left such a number of features. If it is the float value from 0 to 1, the features will be selected automatically in order to save this amount of information. Print the size of the resulting matrix
End of explanation
"""


# %%time


# X_train, X_test, y_train, y_test = train_test_split(df_new, Y,
#                                                     test_size=0.3,
#                                                     random_state=42)

# scaled_data = StandardScaler().fit_transform(X_train)
# scaled_X = pd.DataFrame(scaled_data, index=X_train.index, columns=X_train.columns)

# clf_logreg = LogisticRegression(random_state=42, C=0.05, penalty='l2')

# parameters_grid = {'C': [0.05]}
# columns = X_train.columns
# important_features = []
# features_scores = []


# np.random.seed(123)


# for j in range(len(columns)-1):
#     print('{}\{}'.format(j, len(columns)))
#     col_for_del = []
#     scores = []
#     for i in columns:
#         cols = columns[columns != i]
#         clf = GridSearchCV(clf_logreg, cv=3,
#                            param_grid=parameters_grid,
#                            scoring='f1_macro',
#                            verbose=0,
#                            n_jobs=6).fit(scaled_X[cols], y_train)

#         scores.append(clf.best_score_)
#     max_col = columns[np.argmax(scores)]
#     print(max_col, clf.best_score_)
#     important_features.append(max_col)
#     features_scores.append(max(scores))
#     columns = columns[columns != max_col]
#     print(j, '\r', end='')
# features_scores.append(0)
# important_features.append(columns[0])

# feature_importances = pd.DataFrame({'features': important_features, 
#                                     'feature_importances': features_scores,
#                                     'iteration': range(len(features_scores))})

# feature_importances = feature_importances.sort_values('iteration', ascending=True)
# feature_importances.to_csv(FEATURE_IMPORTANCE_F, index=False)

feature_importances = pd.read_csv(FEATURE_IMPORTANCE_F)

plt.figure(figsize=[15, 6])
plt.title('Brute force')
plt.plot(feature_importances.iteration[:-1], 
         feature_importances.feature_importances[:-1], '-o')
plt.ylabel('f1_macro')
plt.xlabel('iteration')
plt.tight_layout()
plt.show()

"""
Explanation: <a id='greedy_selection'></a>
4.2. Greedy selection
Now we will use the most reliable, but very resource-intensive method - iterative brute force feature selection. We will remove one feature at a time, predict the target variable and see how the quality metric changes: the more the metric decreases after removing the feature, the more important this feature is. More than that, sometimes the metric might even improve.
This method takes a lot of time. Therefore, we won't calculate it now, but you can see the code and its results below.
End of explanation
"""


feature_importances.features[::-1][:20]

"""
Explanation: Top 20 of the most significant features:
End of explanation
"""


feature_importances[feature_importances.feature_importances == feature_importances.feature_importances.max()]

max_id = feature_importances[
    feature_importances.feature_importances == feature_importances.feature_importances.max()
].index[0]

print('We can leave {} features'.format(len(feature_importances)-max_id))

selected_features = feature_importances.features[max_id:]

"""
Explanation: Let's find at what stage was the best accuracy.
End of explanation
"""


# Greedy features

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# X_train, X_test, y_train, y_test = train_test_split(df_new[ ? ], Y,                                                   
#                                                     test_size=0.3, 
#                                                     random_state=42)

pipe = Pipeline([('scale', StandardScaler()),
                 ('clf', LogisticRegression(random_state=42))])

params = {
    'clf__C': [0.01, 0.05, 0.1, 0.5, 0.9, 0.99],
    'clf__penalty': ['l1', 'l2', 'elasticnet']
} 

np.random.seed(123)

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# with ? :
#     clf = GridSearchCV( ? ).fit( ? )
    
learning_time = timer.elapsed_time

with timer:
    predict = clf.predict(X_test)
    
predict_time = timer.elapsed_time

calculate_metrics(predict, y_test)

results = results.append({'method': 'Greedy selection',
                          'model': 'LR',
                          'val score': clf.best_score_,
                          'test score': f1_score(clf.predict(X_test), y_test, average='macro'), 
                          'learning time': learning_time,
                          'predict time': predict_time},
                          ignore_index=True)

# PCA features

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# X_train, X_test, y_train, y_test = train_test_split( ? )

np.random.seed(123)

pipe = Pipeline([('scale', StandardScaler()),
                 ('pca', PCA(n_components=0.99)),
                 ('clf', LogisticRegression(random_state=42))])

##### Implement this part of the code #####
raise NotImplementedError("Code not implemented, follow the instructions.")
# with ?
#     clf = ?    

learning_time = timer.elapsed_time

with timer:
    predict = clf.predict(X_test)
    
predict_time = timer.elapsed_time    

calculate_metrics(clf.predict(X_test), y_test)

results = results.append({'method': 'PCA features',
                          'model': 'LR',
                          'val score': clf.best_score_,
                          'test score': f1_score(clf.predict(X_test), y_test, average='macro'), 
                          'learning time': learning_time,
                          'predict time': predict_time},
                          ignore_index=True)

"""
Explanation: <a id='4_final_prediction'></a>
4.3. Final prediction
Task 19 Train the model on the features obtained by the greedy selection method and the PCA method. Add the results to the table.
End of explanation
"""


results

"""
Explanation: <a id='conclusions'></a>
4.4. Conclusion
End of explanation
"""


# from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
# from sklearn.linear_model import LogisticRegression
# from catboost import CatBoostClassifier
# from lightgbm import LGBMClassifier

# log_params = {'clf__C': uniform(0.01, 0.99),
#               'clf__penalty': ['l1', 'l2']} 

# rf_params = {'clf__n_estimators': randint(100, 1000),
#              'clf__max_depth': randint(1,10),
#              'clf__min_samples_leaf': randint(1, 10),
#              'clf__min_samples_split': randint(2, 11)} 

# gb_params = {'clf__learning_rate': uniform(0.01, 0.49),
#              'clf__n_estimators': randint(100, 700),
#              'clf__max_depth': randint(1, 10),
#              'clf__subsample': uniform(0.6, 0.4)}   

# models_dict = {'Logistic Regression':[log_params, LogisticRegression(random_state=42)],
#                'Random Forest':[rf_params, RandomForestClassifier(random_state=42)],
#                'Gradient Boosting': [gb_params, GradientBoostingClassifier(random_state=42)],
#                'XGB':[gb_params, xgb.XGBClassifier(random_state=42)],
#                'CatBoost':[gb_params, CatBoostClassifier(bootstrap_type='Bernoulli', random_state=42, verbose=0)],
#                'LightBoost':[gb_params, LGBMClassifier(random_state=42)]}

# ext_results = pd.DataFrame(columns=['method', 'model', 'val score', 'test score', 'learning time', 'predict time'])

# def custom_pipe(clf, pca='False'):
#     if pca:
#         return  Pipeline([('clf', clf)])
#     else:
#         return Pipeline([('scale', StandardScaler()),
#                          ('clf', clf)])

# X_train_gs, X_test_gs = train_test_split(df_new[selected_features],test_size=0.3, random_state=42)
# X_train_all, X_test_all = train_test_split(df_new,test_size=0.3, random_state=42)

# pca_model = PCA(n_components=0.99)
# scaler = StandardScaler()
# X_train_pca = pca_model.fit_transform(scaler.fit_transform(X_train_all))
# X_test_pca = pca_model.transform(scaler.transform(X_test_all))

# data_dict = {'Manual features':(X_train_all, X_test_all),  
#              'Greedy selection':(X_train_gs, X_test_gs),
#              'PCA':(X_train_pca, X_test_pca)}

# np.random.seed(123)
# for method, (X_train, X_test) in data_dict.items():
#     print(f'Method: {method}\n')
    
#     if method == "PCA":
#         pca_flag = True
#     else:
#         pca_flag = False
    
#     for clf_name, vals in models_dict.items():
#         print(clf_name, '\n')

#         pipe = custom_pipe(vals[1], pca_flag)
        
#         with timer:
#             clf = RandomizedSearchCV(pipe, cv=3, random_state=123, 
#                                      param_distributions=vals[0], 
#                                      n_jobs=4,
#                                      verbose=1,
#                                      n_iter=100,
#                                      scoring='f1_macro').fit(X_train, y_train)

#         print(clf.best_params_ , '\n')

#         learning_time = int(timer.elapsed_time / 300)

#         with timer:
#             predict = clf.predict(X_test)

#         predict_time = timer.elapsed_time

#         ext_results = ext_results.append({'method': method,
#                                           'model': clf_name,
#                                           'val score': clf.best_score_,
#                                           'test score': f1_score(predict, y_test, average='macro'), 
#                                           'learning time': learning_time,
#                                           'predict time': predict_time},
#                                           ignore_index=True)

# ext_results.to_csv(BONUS_F, index=False)

ext_results = pd.read_csv(BONUS_F)

ext_results

"""
Explanation: The best result on the validation set was shown by a model trained on the features obtained by the greedy selection method. Because some features can decrease the accuracy, removing them allows you to get a better result than on the original dataset. However, on the test set, the greedy algorithm showed slightly worse results than the full dataset. This could happen due to the fact that the optimal set of features was selected for the validation set which in our case might have a slightly different distribution and the dataset isn't particularly large. 
The model trained on the features obtained with the PCA method provides an accuracy comparable to the results provided by the models trained on the raw data, but at the same time it learns much faster and the difference between val and test scores is negligible.
We can't say that some approach is definitely better than others, it depends on the task. In some cases, greedy selection can significantly improve the result, but the time it takes can be extremely long and it can overfit. Using the principal component analysis allows you quickly and easily reduce the dimensionality, which can drastically speed up the training and prediction processes, however there may be a loss in accuracy. There are a lot of other algorithms to perform these actions, from this point it's for you to decide which result is more suitable for you.
<a id='bonus'></a>
5. Bonus
Below you can find an example of comparative training of different models on different sets of features with a wide random selection of hyperparameters. Training takes a lot of time, so we have commented the code. You can load and view the results at the very end of the script. These results were obtained on a slightly different set of features, so they may differ from the results obtained in the workshop.
End of explanation
"""


