<jupyter_start><jupyter_text>Krafthack 7-8. february 2022<jupyter_code>%load_ext autoreload
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Lasso
from sklearn.metrics import mean_absolute_percentage_error as mape
from tqdm.notebook import tqdm

from xgboost import XGBRegressor
from catboost import CatBoostRegressor

import tensorflow as tf
from keras import optimizers, Sequential
from keras.models import Model
from keras.layers import Dense, Input, Activation

import matplotlib.pyplot as plt

from utils.preprocessing import get_timeslots, get_temporal_lookback_features, get_temporal_lookback_df, add_hour_feature, add_seconds_operational<jupyter_output><empty_output><jupyter_text>Import data<jupyter_code>df_train = pd.read_parquet('data/input_dataset-2.parquet')
df_test = pd.read_parquet('data/prediction_input.parquet')<jupyter_output><empty_output><jupyter_text>Clean dataset<jupyter_code># Extract relevant features
cols_keep = list(df_test.columns) + [f'Bolt_{i}_Tensile' for i in range(1,7)]
df_train = df_train[cols_keep]

# Remove rows that contain any missing values
df_train = df_train.dropna()
# Combine both datasets before doing feature engineering
df_full = pd.concat([df_train, df_test], axis=0)
df_full.tail().T<jupyter_output><empty_output><jupyter_text>Feature Engineering- Log-transform- Signal-analysis (derivatives, Fourier transform, power, etc)- Temporal features (day, month, holiday, etc)- Sequencing- Onehhot encoding of categorical features  Add temporal features<jupyter_code>df_full = add_hour_feature(df_full)
df_full = add_seconds_operational(df_full)
df_full['time_weekday'] = df_full.index.dayofweek
df_full.tail().T<jupyter_output><empty_output><jupyter_text>Handle categorical feature<jupyter_code>def get_mode_as_dummy(df):
    # Make "mode" into dummy variable
    y = pd.get_dummies(df["mode"], prefix="Mode")
    df = df.join(y)
    df.drop("mode", inplace=True, axis=1)
    return df

df_full = get_mode_as_dummy(df_full)
df_full.tail()<jupyter_output><empty_output><jupyter_text>Get look-back features<jupyter_code># TODO: Maybe change this for some aggregated features instead
# columns = [
#     'Unit_4_Power',
#     'Unit_4_Reactive Power',
#     'Turbine_Guide Vane Opening',
#     'Turbine_Pressure Drafttube',
#     'Turbine_Pressure Spiral Casing',
#     'Turbine_Rotational Speed'
#     ]


# df_timeslots_list = get_timeslots(df_full)
# df_full_with_lookback = get_temporal_lookback_df(df_timeslots_list, cols=columns, window_size=30, steps=5)<jupyter_output><empty_output><jupyter_text>Split data into train-validate-test<jupyter_code># df_train_new = df_full_with_lookback[df_train.index[0]:df_train.index[-1]].dropna()
# df_test_new  = df_full_with_lookback[df_test.index[0]:df_test.index[-1]]

df_train_new = df_full[df_train.index[0]:df_train.index[-1]].dropna()
df_test_new  = df_full[df_test.index[0]:df_test.index[-1]]
df_test_new['1971-01-31 06:47:00':].head(5)
df_train_new.tail(2)
df_test_new.head(2)
print(f"df_train: {df_train.shape}")
print(f"df_train_new: {df_train_new.shape}")
print(f"df_test: {df_test.shape}")
print(f"df_test_new: {df_test_new.shape}")
PCT_SPLIT = 0.7
labels = [f"Bolt_{i}_Tensile" for i in range(1,7)]

X_train_full = df_train_new.drop(labels, axis=1)  # Official test set
y_train_full = df_train_new[labels]  # Official test labels

X_test_full = df_test_new.drop(labels, axis=1)  # Official test set

X_train_train = X_train_full[:int(PCT_SPLIT*len(X_train_full))]  # Private test set
y_train_train = y_train_full[:int(PCT_SPLIT*len(y_train_full))]  # Private test labels
X_train_val = X_train_full[int(PCT_SPLIT*len(X_train_full)):]  # Private validation set
y_train_val = y_train_full[int(PCT_SPLIT*len(y_train_full)):]  # Private validation labels

print(f"X_train_full:   {X_train_full.shape}")
print(f"y_train_full:   {y_train_full.shape}")
print(f"X_test_full:    {X_test_full.shape}")

print()

print(f"X_train_train:  {X_train_train.shape}")
print(f"y_train_train:  {y_train_train.shape}")

print()

print(f"y_train_val:    {y_train_val.shape}")
print(f"X_train_val:    {X_train_val.shape}")<jupyter_output>X_train_full:   (1750000, 11)
y_train_full:   (1750000, 6)
X_test_full:    (226364, 11)

X_train_train:  (1225000, 11)
y_train_train:  (1225000, 6)

y_train_val:    (525000, 6)
X_train_val:    (525000, 11)<jupyter_text>Save datasets<jupyter_code>X_train_full.to_pickle('data/X_train_full.pkl')
y_train_full.to_pickle('data/y_train_full.pkl')

X_test_full.to_pickle('data/X_test_full.pkl')

X_train_train.to_pickle('data/X_train_train.pkl')
y_train_train.to_pickle('data/y_train_train.pkl')

X_train_val.to_pickle('data/X_train_val.pkl')
y_train_val.to_pickle('data/y_train_val.pkl')<jupyter_output><empty_output><jupyter_text>Load datasets<jupyter_code>del df_full, df_test, df_train, df_test_new, df_train_new
X_train_full = pd.read_pickle('data/X_train_full.pkl')
y_train_full = pd.read_pickle('data/y_train_full.pkl')

X_test_full = pd.read_pickle('data/X_test_full.pkl')

X_train_train = pd.read_pickle('data/X_train_train.pkl')
y_train_train = pd.read_pickle('data/y_train_train.pkl')

X_train_val = pd.read_pickle('data/X_train_val.pkl')
y_train_val = pd.read_pickle('data/y_train_val.pkl')

labels = [f"Bolt_{i}_Tensile" for i in range(1,7)]<jupyter_output><empty_output><jupyter_text>Scaling<jupyter_code>scaler_X = StandardScaler()
scaler_y = StandardScaler()

X_train_scaled = pd.DataFrame(
    scaler_X.fit_transform(X_train_train),
    index = X_train_train.index,
    columns = X_train_train.columns
    )
X_train_val_scaled = pd.DataFrame(
    scaler_X.transform(X_train_val),
    index = X_train_val.index,
    columns = X_train_val.columns
    )
y_train_scaled = pd.DataFrame(
    scaler_y.fit_transform(y_train_train),
    index = y_train_train.index,
    columns = y_train_train.columns
    )
y_train_val_scaled = pd.DataFrame(
    scaler_y.transform(y_train_val),
    index = y_train_val.index,
    columns = y_train_val.columns
    )<jupyter_output><empty_output><jupyter_text>Train Model<jupyter_code>models = {}
model_types = {}
model_configs = {}
model_fit_config = {}

def train_model(models, model_type, model_config, model_fit_config, model_name, labels, X, y):
    models[model_name] = {label: model_type(**model_config) for label in labels}

    for label in tqdm(labels):
        models[model_name][label].fit(X, y[label], **model_fit_config)<jupyter_output><empty_output><jupyter_text>Linear Regression<jupyter_code>MODEL_NAME = 'linreg'
model_types[MODEL_NAME] = LinearRegression
model_configs[MODEL_NAME] = {}
model_fit_config[MODEL_NAME] = {}

train_model(
    models,
    model_type=model_types[MODEL_NAME],
    model_config=model_configs[MODEL_NAME],
    model_fit_config=model_fit_config[MODEL_NAME],
    model_name=MODEL_NAME,
    labels=labels,
    X=X_train_scaled,
    y=y_train_scaled
)<jupyter_output><empty_output><jupyter_text>XGBoost<jupyter_code>MODEL_NAME = 'xgboost'
model_types[MODEL_NAME]=XGBRegressor
model_configs[MODEL_NAME] = dict(
    booster="gbtree",
    learning_rate=0.2,
    gamma=0.1,
    max_depth=6,
    eval_metric="mae")
model_fit_config[MODEL_NAME] = {}

train_model(
    models,
    model_type=model_types[MODEL_NAME],
    model_config=model_configs[MODEL_NAME],
    model_fit_config=model_fit_config[MODEL_NAME],
    model_name=MODEL_NAME,
    labels=labels,
    X=X_train_scaled,
    y=y_train_scaled
)<jupyter_output><empty_output><jupyter_text>CatBoost<jupyter_code>MODEL_NAME = 'catboost'
model_types[MODEL_NAME]=CatBoostRegressor
model_configs[MODEL_NAME] = dict(iterations=400)
model_fit_config[MODEL_NAME] = dict(verbose=False)


train_model(
    models,
    model_type=model_types[MODEL_NAME],
    model_config=model_configs[MODEL_NAME],
    model_fit_config=model_fit_config[MODEL_NAME],
    model_name=MODEL_NAME,
    labels=labels,
    X=X_train_scaled,
    y=y_train_scaled
)<jupyter_output><empty_output><jupyter_text>Sequential modelsA linear model, followed by a random forest model that is trained on the residuals of the first model<jupyter_code>class SequentialModels():
    def __init__(self, model_seq, fit_params=None):
        ''' Instantiate with a list of models and corresponding parameters '''
        self.model_seq = model_seq
        self.fit_params = fit_params

    def fit(self, X, y):
        ''' Fits a sequential set of models'''
        X_new = X.copy()
        y_new = y.to_numpy()
        for model, fit_param in tqdm(zip(self.model_seq, self.fit_params), total=len(self.model_seq)):
            model.fit(X_new, y_new, **fit_param)
            y_new = y_new - model.predict(X_new)

    def predict(self, X):
        ''' Predicts using a sequential set of models '''
        y_hat = np.zeros(shape=(X.shape[0],))
        for model in tqdm(self.model_seq, total=len(self.model_seq)):
            y_hat = y_hat + model.predict(X)
        return y_hat
MODEL_NAME = 'sequential'
model_types[MODEL_NAME] = SequentialModels
model_configs[MODEL_NAME] = dict(
    model_seq= [
        Lasso(),
        CatBoostRegressor(
            iterations=50
        )
    ],
    fit_params=[
        {},
        {'verbose':False}
    ])
model_fit_config[MODEL_NAME] = {}

train_model(
    models,
    model_type=model_types[MODEL_NAME],
    model_config=model_configs[MODEL_NAME],
    model_fit_config=model_fit_config[MODEL_NAME],
    model_name=MODEL_NAME,
    labels=labels,
    X=X_train_scaled,
    y=y_train_scaled
)<jupyter_output><empty_output><jupyter_text>Cross Validation- [Special methods for time-series data](https://medium.com/@soumyachess1496/cross-validation-in-time-series-566ae4981ce4)  Hyperparameter tuning- [Sklearn](https://scikit-learn.org/stable/modules/grid_search.html)- [Nevergrad](https://facebookresearch.github.io/nevergrad/)- [Keras Tuner](https://www.tensorflow.org/tutorials/keras/keras_tuner)  Predict<jupyter_code># labels = [f"Bolt_{i}_Tensile" for i in range(1,7)]
y_preds_scaled = {}
y_preds = {}


for model_name, model in tqdm(models.items(), total=len(models)):
    # For each model type
    y_preds[model_name] = {}
    for label in tqdm(labels):
        # For each sub-model specialized for a unique label column
        y_preds[model_name][label] = model[label].predict(X_train_val_scaled)

    y_preds_scaled[model_name] = pd.DataFrame(y_preds[model_name])
    y_preds[model_name] = pd.DataFrame(scaler_y.inverse_transform(y_preds_scaled[model_name].to_numpy()),
                                       index=y_preds_scaled[model_name].index,
                                       columns=y_preds_scaled[model_name].columns
                                       )<jupyter_output><empty_output><jupyter_text>Visualize results  Plot predicted and true values<jupyter_code>fig, axs = plt.subplots(6, len(y_preds), sharey=True, sharex=True, figsize=(20,10))
fig.suptitle('Predicted and true values', fontsize=16)
for j, (model_name, y_hats) in tqdm(enumerate(y_preds.items()), total=len(y_preds)):
    for i, col in tqdm(enumerate(y_hats.columns), total=y_hats.shape[1]):
        axs[i,j].plot(y_train_val[col].to_numpy(), label='true')
        axs[i,j].plot(y_hats[col].to_numpy(), label='predicted')
        axs[i,j].set_title(f"{model_name} {col}")<jupyter_output><empty_output><jupyter_text>Plot residuals<jupyter_code>fig, axs = plt.subplots(6, len(y_preds), sharey=True, sharex=True, figsize=(20,10))
fig.suptitle('Residuals', fontsize=16)
for j, (model_name, y_hats) in tqdm(enumerate(y_preds.items()), total=len(y_preds)):
    for i, col in tqdm(enumerate(y_hats.columns), total=y_hats.shape[1]):
        axs[i,j].plot(y_hats[col] - y_train_val[col].to_numpy())
        axs[i,j].set_title(f"{model_name} {col}")<jupyter_output><empty_output><jupyter_text>Scoring- Good metrics for temporal data?- Depends on competition metric<jupyter_code>score = {}

for model_name, model in y_preds.items():
    score[model_name] = {}
    
    for label in labels:
        score[model_name][label] = mape(y_train_val[label], y_preds[model_name][label])

scores = pd.DataFrame(score).T
scores.rename({x: f'MAPE {x}' for x in scores.columns}, axis=1)
scores['Avg MAPE'] = scores.mean(axis=1)
scores = scores.sort_values(by='Avg MAPE', ascending=True)
BEST_MODEL = scores.index[0]

scores
scores.plot(kind='bar', figsize=(15,10), title="MAPE scores for each model and each label")<jupyter_output><empty_output><jupyter_text>Model explanation- Explainable model ([interpretml](https://github.com/interpretml/interpret))- Certainty score- [LIME](https://github.com/marcotcr/lime)- [SHAP](https://github.com/slundberg/shap)  Feature Importance  Linear Regression<jupyter_code>pd.DataFrame({label: models['linreg'][label].coef_ for label in labels}, index=X_test_full.columns).plot(kind='bar', figsize=(20,5), title='Coefficients for the linear regression model')<jupyter_output><empty_output><jupyter_text>CatBoost<jupyter_code>pd.DataFrame({label: models['catboost'][label].get_feature_importance() for label in labels}, index=X_test_full.columns).plot(kind='bar', figsize=(20,5), title="Feature importance for the CatBoost model")
import shap
shap.initjs()

explainer = shap.TreeExplainer(models['catboost'][labels[0]])
shap_values = explainer.shap_values(X_test_full_scaled)

# visualize the first prediction's explanation
shap.force_plot(explainer.expected_value, shap_values[0,:],X_test_full_scaled[0,:])<jupyter_output><empty_output><jupyter_text>The above explanation shows features each contributing to push the model output from the base value (the average model output over the training dataset we passed) to the model output. Features pushing the prediction higher are shown in red, those pushing the prediction lower are in blue.  Prepare competition submission  Train model on whole training set<jupyter_code># Scale full training set
scaler_X_full = StandardScaler()
scaler_y_full = StandardScaler()

X_train_full_scaled = scaler_X.fit_transform(X_train_full)
X_test_full_scaled = scaler_X.transform(X_test_full)
y_train_full_scaled = scaler_y_full.fit_transform(y_train_full)

# Train on full training set
train_model(
    models,
    model_type=model_types[BEST_MODEL],
    model_config=model_configs[BEST_MODEL],
    model_fit_config=model_fit_config[BEST_MODEL],
    model_name=BEST_MODEL,
    labels=labels,
    X=X_train_scaled,
    y=y_train_scaled
)<jupyter_output><empty_output><jupyter_text>Make predictions<jupyter_code>y_final_preds_scaled = {}
y_final_preds = {}

for label in tqdm(labels):
    # For each sub-model specialized for a unique label column
    y_final_preds_scaled[label] = models[BEST_MODEL][label].predict(X_test_full_scaled)

y_final_preds_scaled = pd.DataFrame(y_final_preds_scaled)
y_final_preds = pd.DataFrame(scaler_y.inverse_transform(y_final_preds_scaled.to_numpy()),
                                    index=y_final_preds_scaled.index,
                                    columns=y_final_preds_scaled.columns
                                    )<jupyter_output><empty_output><jupyter_text>Export CSV<jupyter_code>y_final_preds.index = X_test_full.index
y_final_preds.to_csv('submission.csv')

y_final_preds.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11755043038307379, 'fraction_numerical': 0.01071029634010119, 'mean_word_length': 5.7196467991169975, 'pattern_counts': {'":': 0, '<': 96, '<?xml version=': 0, '>': 96, 'https://': 7, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '30012', 'license': 'MIT', 'n_tokens_mistral': 5835, 'n_tokens_neox': 5700, 'n_words': 880, 'path': 'main_notebook.ipynb', 'repo_name': 'viggotw/SopraSteria_Team3', 'size': '15219'}