<jupyter_start><jupyter_text>Load Libraries & Datasets<jupyter_code>import os, io, gc
import numpy as np
import pandas as pd
import random

from scipy.fft import fft
from scipy.signal import hilbert, blackman
from sklearn.metrics import mean_absolute_error
from sklearn.preprocessing import RobustScaler
from sklearn.model_selection import KFold

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras.layers import Input, Dense, Dropout, Activation
from tensorflow.keras.layers import Add, concatenate
from tensorflow.keras.layers import Bidirectional, LSTM
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.utils import plot_model
from tensorflow.keras.models import Model, Sequential, load_model
from tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping, ModelCheckpoint

import matplotlib.pyplot as plt
pd.set_option('display.max_columns',None)

import warnings
warnings.filterwarnings('ignore')
gc.enable()
SEED = 42
def seed_everything(seed=SEED):
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    random.seed(seed)
    tf.random.set_seed(seed)
seed_everything()
train_df = pd.read_csv('../input/ventilator-pressure-prediction/train.csv')
test_df = pd.read_csv('../input/ventilator-pressure-prediction/test.csv')
submission_df = pd.read_csv('../input/ventilator-pressure-prediction/sample_submission.csv')
pressure_unique = np.sort(train_df['pressure'].unique())
len_pressure = len(pressure_unique)
PRESSURE_MIN = pressure_unique[0].item()
PRESSURE_MAX = pressure_unique[-1].item()
PRESSURE_STEP = (pressure_unique[1] - pressure_unique[0]).item()<jupyter_output><empty_output><jupyter_text>Data Processing<jupyter_code>def feature_processing(df, backward_sequence=False):
    
    feature_list_1 = ['time_step', 'u_in', 'u_out']
    feature_list_2 = ['u_in', 'u_out']
    
    ###########################################################################

    # backward sequence
    if backward_sequence:
        df['rank'] = df.groupby(['breath_id'])['time_step'].rank()
        df['neg_rank'] = -1 * df['rank']
        df = df.sort_values(by=['breath_id','neg_rank']).reset_index(drop=True)
        del df['rank'],df['neg_rank']
        gc.collect()
    
    ###########################################################################

    df['one'] = 1
    df['count'] = (df['one']).groupby(df['breath_id']).cumsum()
    
    df['u_in_cumsum'] = df.groupby(['breath_id'])['u_in'].cumsum()
    df['u_in_cummean'] = df['u_in_cumsum'] / df['count']
    
    del df['one'], df['count']
    gc.collect()
    
    ###########################################################################
    
    # cross, delta & area

    df['log_u_in'] = np.log1p(df['u_in'])
    df["cross_u_in"] = df["u_in"] * (1 - df["u_out"])
    df["cross_time_step"] = df["time_step"] * (1 - df["u_out"])
    
    df['cross_1']= df['u_in'] * df['u_out']
    df['cross_2']= df['time_step'] * df['u_out']
    
    df["area_frac"] = df["u_in"] * df["time_step"]
    df["cross_area_frac"] = df["area_frac"] * (1 - df["u_out"])
    df["area_frac_cumsum"] = df.groupby(['breath_id'])["area_frac"].cumsum()
    
    df['time_gap'] = df['time_step'] - df.shift(1).fillna(0)['time_step']
    df['u_in_gap'] = df['u_in'] - df.shift(1).fillna(0)['u_in']
    df['u_in_rate'] = df['u_in_gap'] / df['time_gap']
    
    df.loc[list(range(0, len(df), 80)), 'time_gap'] = 0
    df.loc[list(range(0, len(df), 80)), 'u_in_gap'] = 0
    df.loc[list(range(0, len(df), 80)), 'u_in_rate'] = 0
    
    df['area_1'] = df['u_in'] * df['time_step']
    df['area_2'] = df['u_in'] * df['time_gap']
    df['area_timestep_cumsum'] = df.groupby(['breath_id'])['area_1'].cumsum()
        
    df['air_flow_rate'] = df['u_out'] - (df['u_in']/100)
    df['air_flow_area'] = df['air_flow_rate'] * df['time_step']
    
    df['time_step_diff_1'] = df.groupby(['breath_id'])['time_step'].diff(1).fillna(0)
    df['time_step_diff_1_r'] = df.groupby(['breath_id'])['time_step'].diff(-1).fillna(0)
    
    df['delta_1'] = df['time_step_diff_1'] * df['u_in']
    df['delta_2'] = df['time_step_diff_1_r'] * df['u_in']
    
    df['area_1'] = df.groupby(['breath_id'])['delta_1'].cumsum()
    df['area_2'] = df.groupby(['breath_id'])['delta_2'].cumsum()
    df['area_delta_cumsum'] = df.groupby(['breath_id'])['area_1'].cumsum()
    
    df['max_to_cumsum_u_in_per_breath_id'] = df.groupby(['breath_id'])['u_in'].transform('max') - df['u_in_cumsum']
    
    ###########################################################################
    
    # vf: approximation for rate of change in volume at a particular time stamp
    # vt: approximation for total lungs volume at a particular time stamp
    # source: https://www.kaggle.com/c/ventilator-pressure-prediction/discussion/281299
    
    df['vt'] = 0
    df['exponent'] = (-df['time_step']) / (df['R'] * df['C']) 
    df['factor'] = np.exp(df['exponent'])
    df['v1'] = (df['u_in'] * df['R']) / df['factor']
    df['vf'] = (df['u_in_cumsum'] * df['R']) / df['factor']
    df.loc[df['time_step'] != 0, 'vt'] = df['area_timestep_cumsum']/(df['C'] * (1 - df['factor']))
    df['v'] = df['vf'] + df['vt']
        
    ###########################################################################
    
    # lags, difference, and rolling
    
    lags = 3
    for lag in range(1, lags+1):
        for feature in feature_list_1:
            ## lag 
            df[f'{feature}_lag_{lag}'] = df.groupby(['breath_id'])[feature].shift(lag).fillna(0)
            ## inverse lag
            df[f'{feature}_lag_inverse_{lag}'] = df.groupby(['breath_id'])[feature].shift(-lag).fillna(0)
            ## diff lag
            # df[f'{feature}_lag_diff_{lag}'] = df[feature] - df[f'{feature}_lag_{lag}']
            # df[f'{feature}_lag_diff_{lag}'] = df[f'{feature}_lag_diff_{lag}'].fillna(0)
            ## diff inverse lag
            # df[f'{feature}_lag_inverse_diff_{lag}'] = df[feature] - df[f'{feature}_lag_inverse_{lag}']
    
    diff = 3
    for diff in range(1, diff+1):
        df[f'u_in_diff_{diff}'] = df.groupby(['breath_id'])['u_in'].diff(diff).fillna(0)

    # df['u_in_diff_1'] = df.groupby(['breath_id'])['u_in'].diff(1).fillna(0)
    
    """
    lags = 3
    for lag in range(1, lags+1):
        for feature in feature_list_2:
            # breath_id lag 
            df[f'breath_id_lag_{lag}'] = df['breath_id'].shift(lag).fillna(0)
            # breath_id same lag
            df[f'breath_id_lag_{lag}_same'] = np.select([df[f'breath_id_lag_{lag}'] == df['breath_id']],[1],0)
            # breath_id and feature_list_2
            df[f'breath_id_{feature}_lag_{lag}'] = df[feature].shift(lag).fillna(0)
            df[f'breath_id_{feature}_lag_{lag}'] = df[f'breath_id_{feature}_lag_{lag}'] * df[f'breath_id_lag_{lag}_same']
            del df[f'breath_id_lag_{lag}_same'], df[f'breath_id_{feature}_lag_{lag}']
    """
    df['mean_u_out_per_breath_id'] = df.groupby(['breath_id'])['u_out'].transform('mean')
    df['breath_id_u_in_max'] = df.groupby(['breath_id'])['u_in'].transform('max')
    df['breath_id_u_in_diff_max'] = df.groupby(['breath_id'])['u_in'].transform('max') - df['u_in']
    df['breath_id_u_in_diff_max'] = df.groupby(['breath_id'])['u_in'].transform('mean') - df['u_in']
    
    """ 
    windows = [8, 16, 32]
    for feature in feature_list_:
        for window in windows:
            df[f'{feature}_rolling_mean_{window}'] = df.groupby('breath_id')[feature].rolling(window).mean().reset_index(drop=True)
            df[f'{feature}_rolling_min_{window}'] = df.groupby('breath_id')[feature].rolling(window).min().reset_index(drop=True)
            df[f'{feature}_rolling_max_{window}'] = df.groupby('breath_id')[feature].rolling(window).max().reset_index(drop=True)
            df[f'{feature}_rolling_std_{window}'] = df.groupby('breath_id')[feature].rolling(window).std().reset_index(drop=True)
            df[f'{feature}_rolling_std_{window}'] = df.groupby('breath_id')[feature].rolling(window).sum().reset_index(drop=True)
    """
        
    ###########################################################################
    
    # Features based on aggregations over R, C, rank and rounded u_in value (f1 - f6)
    # Source: https://www.kaggle.com/l0glikelihood/0-1093-single-public-lb
    
    df['sum_per_breath'] = df.groupby(['breath_id'])['u_in'].transform('sum')
    df['rounded_u_in'] = df['u_in'].round(0)
    df['rank'] = df.groupby(['breath_id'])['time_step'].rank()
    df['uid'] = df['R'].astype(str)+'_' + df['C'].astype(str) + '_' + df['rounded_u_in'].astype(str) + '_' + df['rank'].astype(str)

    # max, min, mean, count values of u_in for each uid
    df['uid_count'] = df.groupby(['uid'])['uid'].transform('count') 
    df['f1'] = df.groupby(['uid'])['u_in'].transform('mean')
    df['f2'] = df.groupby(['uid'])['u_in'].transform('min')
    df['f3'] = df.groupby(['uid'])['u_in'].transform('max')

    # difference between the current value of u_in and its mean, min and max values within the uid
    df['f4'] = df['u_in'] - df.groupby(['uid'])['u_in'].transform('mean')
    df['f5'] = df['u_in'] - df.groupby(['uid'])['u_in'].transform('min')
    df['f6'] = df['u_in'] - df.groupby(['uid'])['u_in'].transform('max')
    
    del df['rounded_u_in'],df['rank'],df['uid']
        
    ###########################################################################
    
    """
    # spectral features
    # source: https://www.kaggle.com/lucasmorin/spectral-analysis-feature-engineering
    
    ffta = lambda x: np.abs(fft(np.append(x.values,x.values[0]))[:80])
    ffta.__name__ = 'ffta'

    fftw = lambda x: np.abs(fft(np.append(x.values,x.values[0])*w)[:80])
    fftw.__name__ = 'fftw'

    N = 80
    w = blackman(N+1)

    df['fft_u_in'] = df.groupby('breath_id')['u_in'].transform(ffta)
    df['fft_u_in_w'] = df.groupby('breath_id')['u_in'].transform(fftw)
    df['analytical'] = df.groupby('breath_id')['u_in'].transform(hilbert)
    df['envelope'] = np.abs(df['analytical'])
    df['phase'] = np.angle(df['analytical'])
    df['unwrapped_phase'] = df.groupby('breath_id')['phase'].transform(np.unwrap)
    df['phase_shift1'] = df.groupby('breath_id')['unwrapped_phase'].shift(1).astype(np.float32)
    df['IF'] = df['unwrapped_phase'] - df['phase_shift1'].astype(np.float32)
    df = df.fillna(0)
    
    del df['analytical']
    """
    
    ###########################################################################
    
    # R and C features
    
    df['R_u_in'] = df['u_in'] * df['R']
    df['C_u_in'] = df['u_in'] * df['C']
    
    # mean u_in per R, C, and u_out
    df['mean_u_in_per_R_C_u_out'] = df.groupby(['R','C','u_out'])['u_in'].transform('mean')
    df['diff_mean_u_in_per_R_C_u_out'] = df['u_in'] - df['mean_u_in_per_R_C_u_out']
    df['to_mean_u_in_per_R_C_u_out'] = df.groupby(['breath_id'])['u_in'].transform('mean') - df['mean_u_in_per_R_C_u_out']
    
    # max value of u_in grouped by R, C, and u_out
    df['max_u_in_per_R_C_u_out'] = df.groupby(['R','C','u_out'])['u_in'].transform('max')
    df['diff_max_u_in_per_R_C_u_out'] = df['u_in'] - df['max_u_in_per_R_C_u_out']
    df['to_max_u_in_per_R_C_u_out'] = df.groupby(['breath_id'])['u_in'].transform('max') - df['max_u_in_per_R_C_u_out']
    
    # OHE
    df['R'] = df['R'].astype(str)
    df['C'] = df['C'].astype(str)
    df['R_C'] = df['R'].astype(str) + '_' + df['C'].astype(str)
    df = pd.get_dummies(df)
    
    ###########################################################################
        
    return df
df = pd.concat([train_df,test_df],axis=0,copy=False).reset_index(drop=True)
df = feature_processing(df)
    
train = df.iloc[:len(train_df),:]
test = df.iloc[len(train_df):,:].reset_index(drop=True)
    
del df, train_df, test_df
del test['pressure']
gc.collect()
train_df = train.copy()
targets = train[['pressure']].to_numpy().reshape(-1, 80)
train.drop(['pressure','id', 'breath_id'], axis=1, inplace=True)
u_outs = train[['u_out']].to_numpy().reshape(-1, 80)
test = test.drop(['id', 'breath_id'], axis=1)
# Scaler
RS = RobustScaler(quantile_range=(20.0, 80.0))
RS.fit(train[train['u_out']==0])

train = RS.fit_transform(train)
test = RS.transform(test)

train = train.reshape(-1, 80, train.shape[-1])
test = test.reshape(-1, 80, train.shape[-1])
print(f"train: {train.shape} \ntest: {test.shape} \ntargets: {targets.shape} \nu_outs: {u_outs.shape}")
gc.collect()<jupyter_output><empty_output><jupyter_text>ResBiLSTM Model<jupyter_code># Accelerator Configuration
try:
    tpu = tf.distribute.cluster_resolver.TPUClusterResolver()
    tf.config.experimental_connect_to_cluster(tpu)
    tf.tpu.experimental.initialize_tpu_system(tpu)
    tpu_strategy = tf.distribute.experimental.TPUStrategy(tpu)
    BATCH_SIZE = tpu_strategy.num_replicas_in_sync * 64
    print("Running on TPU:", tpu.master())
    print(f"Batch Size: {BATCH_SIZE}")
    
except ValueError:
    strategy = tf.distribute.get_strategy()
    BATCH_SIZE = 512
    print(f"Running on {strategy.num_replicas_in_sync} replicas")
    print(f"Batch Size: {BATCH_SIZE}")
# Model Configuration
class CFG:
    seed = 42
    VERBOSE = 1
    random_state = 42
    N_FOLDS = 5
    EPOCHS = 200
    BATCH_SIZE = BATCH_SIZE
    factor = 0.5
    patience_1 = 5
    patience_2 = 15
    learning_rate = 1e-3
    weight_decay = 1e-3
    dropout_rate = 0.2
# Custom MAE Loss
def custom_mae_loss(y_true, y_pred, n=80):
    u_out = y_true[:, n: ]
    y = y_true[:, :n ]
    w = 1 - u_out
    mae = w * tf.abs(y - y_pred)
    return tf.reduce_sum(mae, axis=-1) / tf.reduce_sum(w, axis=-1)
def dnn_model():
    x_input = Input(shape=(train.shape[-2:]))
    
    x1 = Bidirectional(LSTM(units=1024, return_sequences=True))(x_input)
    c1 = concatenate([x_input, x1])
    
    x2 = Bidirectional(LSTM(units=512, return_sequences=True))(c1)
    c2 = concatenate([x1, x2])
    
    x3 = Bidirectional(LSTM(units=256, return_sequences=True))(c2)
    c3 = concatenate([x2, x3])
    
    x4 = Bidirectional(LSTM(units=128, return_sequences=True))(c3)
    c4 = concatenate([x3, x4])
    
    x5 = Dense(units=128, activation='selu')(c4)
    x_output = Dense(units=1)(x5)
    
    model = Model(inputs=x_input, outputs=x_output, name='DNN_Model')
    model.compile(optimizer='Adam', loss=custom_mae_loss)
    
    return model
model = dnn_model()
model.summary()
plot_model(model, to_file='dnn_model.png', show_shapes=True, show_layer_names=True)
test_preds = []
history_list = []
oof_true = []
oof_pred = []
with tpu_strategy.scope():
    kf = KFold(n_splits=CFG.N_FOLDS, shuffle=True, random_state=42)
    oof_preds = np.zeros((train.shape[0],train.shape[1]))

    for fold, (train_idx, test_idx) in enumerate(kf.split(train, targets)):
        print('='*25, '>', f'Fold {fold+1}', '<', '='*25)
        
        checkpoint_filepath = f'fold{fold+1}.hdf5'
        X_train, X_valid = train[train_idx], train[test_idx]
        y_train, y_valid = targets[train_idx], targets[test_idx]
        u_out_train, u_out_valid = u_outs[train_idx], u_outs[test_idx] 
            
        lr = ReduceLROnPlateau(monitor="val_loss", factor=CFG.factor, patience=CFG.patience_1,
                               verbose=CFG.VERBOSE)
        es = EarlyStopping(monitor="val_loss", patience=CFG.patience_2, mode="min",
                           restore_best_weights=True, verbose=CFG.VERBOSE)
        sv = ModelCheckpoint(checkpoint_filepath, monitor = 'val_loss', verbose = CFG.VERBOSE,
                             save_best_only = True, save_weights_only = True, mode = 'min')
        
        model = dnn_model()
        history = model.fit(X_train, np.append(y_train, u_out_train, axis =1),
                            validation_data=(X_valid, np.append(y_valid, u_out_valid, axis =1)),
                            epochs=CFG.EPOCHS, batch_size=CFG.BATCH_SIZE, callbacks=[lr,es,sv])        
        history_list += [history]
        
        # predict oof
        y_pred = model.predict(X_valid)
        y_true = y_valid.squeeze().reshape(-1, 1)

        ## inspiratory and expiratory phases
        score = mean_absolute_error(y_true, y_pred.squeeze().reshape(-1, 1))
        print(f'Fold {fold+1} | Overall MAE Score: {score}')
        
        ## inspiratory phase
        oof_true.append(y_true)
        oof_pred.append(y_pred.squeeze().reshape(-1, 1))
        oof_preds[test_idx] = y_pred.reshape(y_pred.shape[0],y_pred.shape[1])
        reshaped_targets = targets.squeeze().reshape(-1,1).squeeze()
        
        score = mean_absolute_error(reshaped_targets,oof_preds.squeeze().reshape(-1,1).squeeze())
        print(f'Fold {fold+1} | Inspiratory MAE Score: {score}')
        
        # predict test
        test_preds.append(model.predict(test).squeeze().reshape(-1, 1).squeeze())
        del X_train, X_valid, y_train, y_valid
        gc.collect()
        
    np.save('test_preds.npy', test_preds)
oof_preds = oof_preds.squeeze().reshape(-1,1).squeeze()
reshaped_targets = targets.squeeze().reshape(-1,1).squeeze()
score = mean_absolute_error(reshaped_targets, oof_preds)
print(f'Overall OOF MAE Score: {score}')
idx = train_df[train_df['u_out']==0].index
train_df['prediction'] = oof_preds
score = mean_absolute_error(train_df.loc[idx,'pressure'],train_df.loc[idx,'prediction'])
print(f'Training Inspiratory MAE Score: {score}')
idx = train_df[train_df['u_out']==1].index
train_df['prediction'] = oof_preds
score = mean_absolute_error(train_df.loc[idx,'pressure'],train_df.loc[idx,'prediction'])
print(f'Training Expiratory MAE Score: {score}')
t = 0
for k in range(CFG.N_FOLDS):
    mae = np.mean(np.abs(oof_pred[k] - oof_true[k]))
    t += mae
    print(f'Fold {k+1} | MAE Validation Score: {mae}')
print(f'Overall CV MAE: {t/CFG.N_FOLDS}')
t = 0
for k in range(CFG.N_FOLDS):
    mae = np.mean(np.abs(oof_preds[k] - oof_true[k]))
    t += mae
    print(f'Fold {k+1} | Inspiratory MAE Score: {mae}')
print(f'Overall Inspiratory MAE Score: {t/CFG.N_FOLDS}')
def plot_hist(hist, with_grid=True):
    plt.figure(figsize=(20,5))
    for i in range(len(hist)):
        plt.plot(hist[i].history["loss"], color='grey')
        plt.plot(hist[i].history["val_loss"], color='green')
    plt.title("")
    plt.ylabel("Mean Absolute Error")
    plt.xlabel("epoch")
    plt.legend(["Training", "Validation"], loc="upper right")
    if with_grid:
        plt.grid(which='major', axis='both')
    plt.show()
plot_hist(history_list)<jupyter_output><empty_output><jupyter_text>Post-Processing<jupyter_code># Ensemble Folds with Mean
submission_df['pressure'] = np.mean(np.vstack(test_preds),axis=0)
submission_df.to_csv('submission_mean.csv', index=False)
# Ensemble Folds with Median
submission_df['pressure'] = np.median(np.vstack(test_preds),axis=0)
submission_df.to_csv('submission_median.csv', index=False)
# Ensemble Folds with Median and Round Prediction
submission_df['pressure'] = np.mean(np.vstack(test_preds),axis=0)
submission_df['pressure'] = np.round((submission_df['pressure'] - PRESSURE_MIN)/PRESSURE_STEP) * PRESSURE_STEP + PRESSURE_MIN
submission_df['pressure'] = np.clip(submission_df['pressure'], PRESSURE_MIN, PRESSURE_MAX)
submission_df.to_csv('submission_median_round.csv', index=False)
# Nearest Neighbor Method
## Mean
submission_df['pressure'] = np.mean(np.vstack(test_preds),axis=0)
submission_df['pressure'] = submission_df['pressure'].map(lambda x: pressure_unique[np.abs(pressure_unique - x).argmin()])
submission_df.to_csv('submission_nn_mean.csv', index=False)

## Median
submission_df['pressure'] = np.median(np.vstack(test_preds),axis=0)
submission_df['pressure'] = submission_df['pressure'].map(lambda x: pressure_unique[np.abs(pressure_unique - x).argmin()])
submission_df.to_csv('submission_nn_median.csv', index=False)
# Mean-Median Method
# Source: https://www.kaggle.com/c/ventilator-pressure-prediction/discussion/282735

def better_than_median(inputs, spread_lim = None, axis=0):
    """
    Compute the mean of the predictions if there are no outliers,
    or the median if there are outliers.

    Parameter: inputs = ndarray of shape (n_samples, n_folds)
    """
    spread = inputs.max(axis=axis) - inputs.min(axis=axis) 
    print(f"Inliers:  {(spread < spread_lim).sum():7} -> compute mean")
    print(f"Outliers: {(spread >= spread_lim).sum():7} -> compute median")
    print(f"Total:    {len(inputs):7}")
    return np.where(spread < spread_lim,
                    np.mean(inputs, axis=axis),
                    np.median(inputs, axis=axis))

submission_df['pressure'] = better_than_median(np.vstack(test_preds), spread_lim = 0.50)
submission_df.to_csv('submission_mixed_50.csv', index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20631456423686892, 'fraction_numerical': 0.013676705950822057, 'mean_word_length': 4.236160487557136, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 21, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '46187', 'license': 'MIT', 'n_tokens_mistral': 7943, 'n_tokens_neox': 7651, 'n_words': 1114, 'path': 'notebook.ipynb', 'repo_name': 'abdelghanibelgaid/nnvpp', 'size': '20619'}