{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# ","metadata":{}},{"cell_type":"markdown","source":"# 1.: Data Preprocessing and Feature Extraction\n## 1.1.: Import Packages and Specify File Paths","metadata":{}},{"cell_type":"code","source":"# Install tsflex and seglearn package (online)\n!pip install tsflex\n!pip install seglearn\n\n# Import packages and functions\nimport numpy as np\nimport pathlib\nimport pandas as pd \nfrom tqdm.auto import tqdm\nfrom sklearn import *\nimport glob\nimport os\nimport gc\nfrom IPython.display import FileLink\nimport sys\nfrom sklearn.model_selection import train_test_split\nfrom seglearn.feature_functions import base_features, emg_features\nfrom tsflex.features import FeatureCollection, MultipleFeatureDescriptors\nfrom tsflex.features.integrations import seglearn_feature_dict_wrapper\n\n# Define important paths\nos.chdir(r'/kaggle/working') \ndataset_path = '/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/'\ndata_path_dataset = '/kaggle/input/fog-features/'\ndata_path_output = '/kaggle/working/'\nmodel_path = '/kaggle/input/ensemble-models/' # '/kaggle/input/regressor-models/' ","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:32:20.094041Z","iopub.execute_input":"2023-06-10T15:32:20.094450Z","iopub.status.idle":"2023-06-10T15:32:47.552447Z","shell.execute_reply.started":"2023-06-10T15:32:20.094420Z","shell.execute_reply":"2023-06-10T15:32:47.551070Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1.2.: Initialize Functions","metadata":{}},{"cell_type":"code","source":"# Reduce Memory Usage\n# Taken and/or adapted from: https://www.kaggle.com/code/arjanso/reducing-dataframe-memory-size-by-65 @ARJANGROEN\n\ndef reduce_memory_usage(df):\n    \"\"\"\n    Reduces the memory usage of a DataFrame by downcasting numeric types and converting object columns to category type.\n    \n    Args:\n        df (pd.DataFrame): The DataFrame to reduce memory usage.\n        \n    Returns:\n        pd.DataFrame: The DataFrame with reduced memory usage.\n    \"\"\"\n    # Get the initial memory usage of the DataFrame\n    start_mem = df.memory_usage().sum() / 1024 ** 2\n    print('Memory usage of the DataFrame is {:.2f} MB'.format(start_mem))\n    \n    # Iterate over each column in the DataFrame\n    for col in df.columns:\n        col_type = df[col].dtype.name\n        \n        # Check if the column type is not datetime or category\n        if col_type not in ['datetime64[ns]', 'category']:\n            \n            # Check if the column type is not object\n            if col_type != 'object':\n                c_min = df[col].min()\n                c_max = df[col].max()\n\n                if str(col_type)[:3] == 'int':\n                    # Downcast integer types based on their minimum and maximum values\n                    if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n                        df[col] = df[col].astype(np.int8)\n                    elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n                        df[col] = df[col].astype(np.int16)\n                    elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n                        df[col] = df[col].astype(np.int32)\n                    elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n                        df[col] = df[col].astype(np.int64)\n\n                else:\n                    # Downcast float types based on their minimum and maximum values\n                    if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n                        df[col] = df[col].astype(np.float16)\n                    elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n                        df[col] = df[col].astype(np.float32)\n            else:\n                # Convert object columns to category type\n                df[col] = df[col].astype('category')\n    \n    # Calculate the final memory usage of the DataFrame\n    mem_usg = df.memory_usage().sum() / 1024 ** 2 \n    print(\"Memory usage reduced to: {:.2f} MB\".format(mem_usg))\n    \n    return df\n\n\n# Read and preprocess data + extract features\n# Adapted from : XXX\ndef reader(f):\n    \"\"\"\n    Reads a CSV file, performs data processing operations, and returns a DataFrame.\n    \n    Args:\n        f (str): File path of the CSV file.\n        tasks (pd.DataFrame): DataFrame containing tasks data.\n        metadata_complex (pd.DataFrame): DataFrame containing complex metadata.\n        fc (FeatureCollection): The feature collection object used for calculating features.\n        \n    Returns:\n        pd.DataFrame: The processed DataFrame.\n    \"\"\"\n    try:\n        # Read the CSV file with the specified index column as 'Time'\n        df = pd.read_csv(f, index_col=\"Time\")\n\n        # Extract the ID and module information from the file path\n        df['Id'] = f.split('/')[-1].split('.')[0]\n        df['Module'] = pathlib.Path(f).parts[-2]\n\n        # Calculate the fractional time index\n        df['Time_frac'] = (df.index / df.index.max()).values\n\n        # Merge with the tasks and metadata_complex DataFrames\n        df = pd.merge(df, tasks[['Id','t_kmeans']], how='left', on='Id').fillna(-1)\n        df = pd.merge(df, metadata_complex[['Id','Subject']+['Visit','Test','Medication','s_kmeans']], how='left', on='Id').fillna(-1)\n\n        # Calculate features using the provided FeatureCollection object (fc)\n        df_feats = fc.calculate(df, return_df=True, include_final_window=True, approve_sparsity=True, window_idx=\"begin\").astype(np.float32)\n\n        # Merge the calculated features with the original DataFrame\n        df = df.merge(df_feats, how=\"left\", left_index=True, right_index=True)\n\n        # Fill missing values with the previous values\n        df.fillna(method=\"ffill\", inplace=True)\n\n        # Modify the ID column to include the index values\n        df['Id'] = df['Id'].astype(str) + '_' + df.index.astype(str)\n\n        return df\n    except:\n        pass","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:32:47.554905Z","iopub.execute_input":"2023-06-10T15:32:47.555860Z","iopub.status.idle":"2023-06-10T15:32:47.580479Z","shell.execute_reply.started":"2023-06-10T15:32:47.555820Z","shell.execute_reply":"2023-06-10T15:32:47.579289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1.3.: Read Auxilary Data","metadata":{}},{"cell_type":"code","source":"# Load the subject and task data from CSV files\nsubjects = pd.read_csv(dataset_path + 'subjects.csv')\ntasks = pd.read_csv(dataset_path + 'tasks.csv')\n\n# Load the sample submission data\nsample_submission = pd.read_csv(dataset_path + 'sample_submission.csv')\n\n# Load the metadata files\ntdcsfog_metadata = pd.read_csv('/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/tdcsfog_metadata.csv')\ndefog_metadata = pd.read_csv('/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/defog_metadata.csv')","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:32:47.582856Z","iopub.execute_input":"2023-06-10T15:32:47.583327Z","iopub.status.idle":"2023-06-10T15:32:47.873020Z","shell.execute_reply.started":"2023-06-10T15:32:47.583283Z","shell.execute_reply":"2023-06-10T15:32:47.871912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1.4.: Preprocess and Store Auxilary Data","metadata":{}},{"cell_type":"code","source":"# Taken and/or adapted from: XXX\n\n# Add a 'Module' column to indicate the source of the metadata\ntdcsfog_metadata['Module'] = 'tdcsfog'\ndefog_metadata['Module'] = 'defog'\n\n# Concatenate the metadata from both sources\nmetadata = pd.concat([tdcsfog_metadata, defog_metadata])\n\n# Calculate the duration of each task\ntasks['Duration'] = tasks['End'] - tasks['Begin']\n\n# Pivot the tasks DataFrame to get the sum of durations for each Id and Task\ntasks = pd.pivot_table(tasks, values=['Duration'], index=['Id'], columns=['Task'], aggfunc='sum', fill_value=0)\n\n# Rename the columns to remove the multi-index and keep only the last character\ntasks.columns = [c[-1] for c in tasks.columns]\n\n# Reset the index of the tasks DataFrame\ntasks = tasks.reset_index()\n\n# Perform K-means clustering on the tasks data using 10 clusters\ntasks['t_kmeans'] = cluster.KMeans(n_clusters=10, random_state=3).fit_predict(tasks[tasks.columns[1:]])\n\n# Fill missing values with 0 and calculate the median for each subject\nsubjects = subjects.fillna(0).groupby('Subject').median()\n\n# Reset the index of the subjects DataFrame\nsubjects = subjects.reset_index()\n\n# Perform K-means clustering on the subjects data using 10 clusters\nsubjects['s_kmeans'] = cluster.KMeans(n_clusters=10, random_state=3).fit_predict(subjects[subjects.columns[1:]])\n\n# Rename the columns of the subjects DataFrame\nsubjects = subjects.rename(columns={\n    'Visit': 's_Visit',\n    'Age': 's_Age',\n    'YearsSinceDx': 's_YearsSinceDx',\n    'UPDRSIII_On': 's_UPDRSIII_On',\n    'UPDRSIII_Off': 's_UPDRSIII_Off',\n    'NFOGQ': 's_NFOGQ'\n})\n\n# Define the list of complex features\ncomplex_featlist = ['Visit', 'Test', 'Medication', 's_Visit', 's_Age', 's_YearsSinceDx', 's_UPDRSIII_On', 's_UPDRSIII_Off', 's_NFOGQ', 's_kmeans']\n\n# Merge the metadata DataFrame with the subjects DataFrame based on the 'Subject' column\nmetadata_complex = metadata.merge(subjects, how='left', on='Subject').copy()\n\n# Convert the 'Medication' column to integer labels using factorize\nmetadata_complex['Medication'] = metadata_complex['Medication'].factorize()[0]\n\n\nsubjects.to_csv('subjects.csv', index=False)\ntasks.to_csv('tasks.csv', index=False)\nmetadata_complex.to_csv('metadata_complex.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:32:47.876327Z","iopub.execute_input":"2023-06-10T15:32:47.876765Z","iopub.status.idle":"2023-06-10T15:32:48.187248Z","shell.execute_reply.started":"2023-06-10T15:32:47.876725Z","shell.execute_reply":"2023-06-10T15:32:48.186421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1.5.: Create Feature Collection","metadata":{}},{"cell_type":"code","source":"# Specify the stride and window size\nwindow_size = 5_000\nstride = 5_000\n\n# Define basic_feats feature descriptors\nbasic_feats = MultipleFeatureDescriptors(\n    functions=seglearn_feature_dict_wrapper(base_features()),\n    series_names=['AccV', 'AccML', 'AccAP'],\n    windows=[window_size],\n    strides=[stride]\n)\n\n# Define emg_feats feature descriptors\nemg_feats = emg_features()\ndel emg_feats['simple square integral']  # Remove 'simple square integral' from emg_feats\n\nemg_feats = MultipleFeatureDescriptors(\n    functions=seglearn_feature_dict_wrapper(emg_feats),\n    series_names=['AccV', 'AccML', 'AccAP'],\n    windows=[window_size],\n    strides=[stride]\n)\n\n# Create a FeatureCollection by combining basic_feats and emg_feats\nfc = FeatureCollection([basic_feats, emg_feats])\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:32:48.189271Z","iopub.execute_input":"2023-06-10T15:32:48.189625Z","iopub.status.idle":"2023-06-10T15:32:48.383257Z","shell.execute_reply.started":"2023-06-10T15:32:48.189594Z","shell.execute_reply":"2023-06-10T15:32:48.382122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1.6.: Create or Load Feature Set","metadata":{}},{"cell_type":"code","source":"# Specify feature extraction characteristics\nextract_features = True # Do you want to extract the features? If features already in a dedicated dataset, they can be loaded in the model section\nvalid = False # Only use valid samples\nreduce = True # Reduce the memory usage of the final feature set\nsafe = False # Safe the final feature set\ndatasets = ['defog', 'tdcsfog'] #, 'notype' #Define what data to use for the feature set\nname = 'train_features' # Define the name of the final dataset\n\n\nif extract_features:\n    train_features = pd.DataFrame()\n    for data in datasets:\n        print(f'Creating {data} features...')\n        \n        # Specify file path and create feature set\n        train_files = glob.glob(dataset_path + 'train/' + data +'/**')\n        features = pd.concat([reader(file) for file in tqdm(train_files)]).fillna(0) #defog_train_files\n        print(f'{data} features created')\n\n        if valid:\n            # Only use valid samples\n            name += '_valid'\n            print('Using only valid samples')\n\n            if 'Valid' in features.columns:\n                features = features[(features['Valid'] == 1) & (features['Task'] == 1)].drop(columns=['Valid', 'Task'])\n        \n        # Drop the 'Valid' and 'Task' column\n        if 'Valid' in features.columns:\n            features = features.drop(columns=['Valid', 'Task'])\n\n        train_features = pd.concat([train_features, features])\n\n        del features\n        gc.collect()\n\n    print('Training features created')\n    train_features.reset_index(drop=True, inplace=True)\n\n    if reduce:\n        # Reduce memory usage of the dataframe\n        print('Reducing memory usage of the DataFrame')\n        name += '_reduced'\n        train_features = reduce_memory_usage(train_features)\n\n    gc.collect()\n    train_features","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:32:48.384728Z","iopub.execute_input":"2023-06-10T15:32:48.385075Z","iopub.status.idle":"2023-06-10T15:46:09.183591Z","shell.execute_reply.started":"2023-06-10T15:32:48.385046Z","shell.execute_reply":"2023-06-10T15:46:09.181960Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1.7.: Split the Feature Set into a Training and Evaluation Set ","metadata":{}},{"cell_type":"code","source":"# Select columns for modeling and prediction excluding unnecessary columns\nfeature_cols = [c for c in train_features.columns if c not in ['Id', 'Subject', 'Module', 'Time', 'StartHesitation', 'Turn', 'Walking', 'Valid', 'Task', 'Event']]\nprediction_cols = ['StartHesitation', 'Turn', 'Walking']\n\n# Split the dataframe into a training and evaluation set\ntrain_set, eval_set = train_test_split(train_features, stratify= train_features[prediction_cols], test_size=0.2, random_state=42)\ndel train_features\ngc.collect()\n\n# Reset the index of the training and evaluation feature sets\ntrain_set.reset_index(drop=True, inplace=True)\neval_set.reset_index(drop=True, inplace=True)\n\nif safe == True:\n    # Store the training and evaluation feature sets\n    print('Safe training features to output folder...')\n    train_set.to_csv('train_set.csv', index=False)\n    print('Training features stored in output')\n    \n    print('Safe evaluation features to output folder...')\n    eval_set.to_csv('eval_set.csv', index=False)\n    print('Evaluation features stored in output')","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:46:09.186435Z","iopub.execute_input":"2023-06-10T15:46:09.186951Z","iopub.status.idle":"2023-06-10T15:50:20.637070Z","shell.execute_reply.started":"2023-06-10T15:46:09.186905Z","shell.execute_reply":"2023-06-10T15:50:20.635723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2.: Train the Ensemble Regressors","metadata":{}},{"cell_type":"markdown","source":"## 2.1.: Import Packages","metadata":{}},{"cell_type":"code","source":"import pickle\nfrom sklearn.model_selection import GroupKFold\nfrom sklearn import metrics\nfrom sklearn.preprocessing import LabelBinarizer\nimport matplotlib.pyplot as plt\nfrom sklearn.multioutput import MultiOutputRegressor\nfrom sklearn.metrics import average_precision_score\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.metrics import roc_curve, auc\nfrom itertools import cycle\n\nfrom xgboost import XGBRegressor\nfrom lightgbm import LGBMRegressor\nfrom sklearn.neural_network import MLPRegressor\nfrom sklearn.gaussian_process import GaussianProcessRegressor\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.ensemble import HistGradientBoostingRegressor\nfrom sklearn.ensemble import AdaBoostRegressor\nfrom sklearn.neighbors import KNeighborsRegressor\nfrom sklearn.ensemble import ExtraTreesRegressor","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:50:20.639016Z","iopub.execute_input":"2023-06-10T15:50:20.639408Z","iopub.status.idle":"2023-06-10T15:50:22.056918Z","shell.execute_reply.started":"2023-06-10T15:50:20.639376Z","shell.execute_reply":"2023-06-10T15:50:22.055883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.2.: Initialize Functions","metadata":{}},{"cell_type":"code","source":"# MultiOutputRegressor\nclass CustomMultiOutputRegressor(MultiOutputRegressor):\n    def fit(self, X, y, eval_set=None, **fit_params):\n        self.estimators_ = [clone(self.estimator) for _ in range(y.shape[1])]\n        \n        for i, estimator in enumerate(self.estimators_):\n            if eval_set:\n                fit_params['eval_set'] = [(eval_set[0], eval_set[1][:, i])]\n            estimator.fit(X, y[:, i], **fit_params)\n        \n        return self\n\n\ndef compute_multiclass_roc_auc_metrics(y_test, y_score):\n\n    n_classes = len(prediction_cols)\n    label_binarizer = LabelBinarizer().fit(train_set[prediction_cols]) # train_set[pcols]\n    y_onehot_test = label_binarizer.transform(y_test)\n    y_onehot_test.shape  # (n_samples, n_classes)\n    \n    # store the fpr, tpr, and roc_auc for all averaging strategies\n    fpr, tpr, roc_auc = dict(), dict(), dict()\n    \n    # Compute micro-average ROC curve and ROC area\n    fpr_micro, tpr_micro, _ = roc_curve(y_onehot_test.ravel(), y_score.ravel())\n    roc_auc_micro = auc(fpr_micro, tpr_micro)\n    \n    \n    for i in range(n_classes):\n        fpr[i], tpr[i], _ = roc_curve(y_onehot_test[:, i], y_score[:, i])\n        roc_auc[i] = auc(fpr[i], tpr[i])\n        \n    fpr_grid = np.linspace(0.0, 1.0, 1000)\n\n    # Interpolate all ROC curves at these points\n    mean_tpr = np.zeros_like(fpr_grid)\n\n    for i in range(n_classes):\n        mean_tpr += np.interp(fpr_grid, fpr[i], tpr[i])  # linear interpolation\n\n    # Average it and compute AUC\n    mean_tpr /= n_classes\n\n    fpr_macro = fpr_grid\n    tpr_macro = mean_tpr\n    roc_auc_macro = auc(fpr_macro, tpr_macro)\n    \n    return fpr, tpr, roc_auc, fpr_micro, fpr_macro, tpr_micro, tpr_macro, roc_auc_micro, roc_auc_macro \n\n\ndef plot_multiclass_ROC(fprs, tprs, fpr_micro_macro, tpr_micro_macro, roc_auc_ovr_micro_macro):\n    \n    micro_macro_names = ['micro', 'macro']\n    micro_macro_colors = ['navy', 'deeppink']\n    class_colors = [\"aqua\", \"darkorange\", \"cornflowerblue\"]\n    n_classes = len(prediction_cols)\n    tprs_micro = []\n    tprs_macro = []\n    fpr_cs = []\n    tpr_cs = []\n    aucs = np.empty([n_classes, N_FOLDS])\n    base_fpr = np.linspace(0, 1, 101)\n\n    plt.figure(figsize=(6, 6))\n    plt.axes().set_aspect('equal', 'datalim')\n\n    # Plot micro / macro ROC\n    for l in range(2):\n        fpr_init = fpr_micro_macro[l]\n        tpr_init = tpr_micro_macro[l]\n        tprs_m = []\n\n        for i in range(N_FOLDS):\n            fpr= fpr_init[i]\n            tpr = tpr_init[i]\n\n            plt.plot(fpr, tpr, micro_macro_colors[l], linestyle='dashed',alpha=0.15)\n            tpr_m = np.interp(base_fpr, fpr, tpr)\n            tpr_m[0] = 0.0\n            tprs_m.append(tpr_m)\n\n        tprs_m = np.array(tprs_m)\n        mean_tprs = tprs_m.mean(axis=0)\n        std = tprs_m.std(axis=0)\n\n        tprs_upper = np.minimum(mean_tprs + std, 1)\n        tprs_lower = mean_tprs - std\n\n        plt.plot(base_fpr, mean_tprs, micro_macro_colors[l],linestyle = 'dashed', label=f\"{micro_macro_names[l]}-average ROC curve (AUC = {np.mean(roc_auc_ovr_micro_macro[l]):.2f} ± {np.std(roc_auc_ovr_micro_macro[l]):.2f})\")\n        plt.fill_between(base_fpr, tprs_lower, tprs_upper, color=micro_macro_colors[l], alpha=0.3)\n\n    # Convert fpr and tpr dict to list\n    for i in range(N_FOLDS):\n        fpr_c = []\n        tpr_c = []\n        fpr_data = list(fprs[i].items())\n        tpr_data = list(tprs[i].items())\n\n        for c in range(3):\n            fpr_c.append(fpr_data[c][1])\n            tpr_c.append(tpr_data[c][1])\n\n        fpr_cs.append(fpr_c)\n        tpr_cs.append(tpr_c)\n\n    # Plot class ROC\n    for c in range(n_classes):\n        tprs_c = []\n\n        for i in range(N_FOLDS): \n            fpr = fpr_cs[i][c]\n            tpr = tpr_cs[i][c]\n\n            plt.plot(fpr, tpr, class_colors[c], alpha=0.15)\n            tpr_c = np.interp(base_fpr, fpr, tpr)\n            tpr_c[0] = 0.0\n            tprs_c.append(tpr_c)\n\n            aucs[c][i]= metrics.auc(fpr, tpr)\n\n        tprs_c = np.array(tprs_c)\n        mean_tprs = tprs_c.mean(axis=0)\n        std = tprs_c.std(axis=0)\n\n        tprs_upper = np.minimum(mean_tprs + std, 1)\n        tprs_lower = mean_tprs - std\n\n        plt.plot(base_fpr, mean_tprs, class_colors[c], label=f\"ROC curve for {prediction_cols[c]} (AUC = {np.mean(aucs[c][:]):.2f} ± {aucs[c][:].std():.2f} )\")\n        plt.fill_between(base_fpr, tprs_lower, tprs_upper, color=class_colors[c], alpha=0.2)\n\n    plt.plot([0, 1], [0, 1],'k:')\n    plt.xlim([-0.01, 1.01])\n    plt.ylim([-0.01, 1.01])\n    plt.ylabel('True Positive Rate')\n    plt.xlabel('False Positive Rate')\n    plt.title(f'{model_description}  \\nOne-vs-Rest multiclass Operating Characteristic')\n    plt.legend()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:50:22.058495Z","iopub.execute_input":"2023-06-10T15:50:22.058822Z","iopub.status.idle":"2023-06-10T15:50:22.092553Z","shell.execute_reply.started":"2023-06-10T15:50:22.058795Z","shell.execute_reply":"2023-06-10T15:50:22.090896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.3.: Setup Regressor Training","metadata":{}},{"cell_type":"code","source":"N_FOLDS = 5\nmodel_type = 'XT'\nValid = False\n\nmatch model_type:\n    case 'LGBM':\n        model = LGBMRegressor()\n    case 'XGB':\n        model = XGBRegressor()\n    case 'ADAB':\n        model = AdaBoostRegressor()\n    case 'HGB':\n        model = HistGradientBoostingRegressor()\n    case 'RF':\n        model = RandomForestRegressor()\n    case 'MLP':\n        model = MLPRegressor()\n    case 'KN':\n        model = KNeighborsRegressor()\n    case 'GP':\n        model = GaussianProcessRegressor()\n    case 'XT':\n        model = ExtraTreesRegressor()","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:50:22.097165Z","iopub.execute_input":"2023-06-10T15:50:22.098127Z","iopub.status.idle":"2023-06-10T15:50:22.109526Z","shell.execute_reply.started":"2023-06-10T15:50:22.098077Z","shell.execute_reply":"2023-06-10T15:50:22.108569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.4.: Load the Training and Evaluation Feature Set if not Created","metadata":{}},{"cell_type":"code","source":"# Load training set\nname = 'train_set'\nif name in globals():\n    pass\nelse:\n    if valid == True:\n        name = name + '_valid' \n        \n    if os.path.isfile(data_path_dataset + name + '.csv'):\n        train_set = pd.read_csv(ata_path_dataset + name + '.csv')\n        train_set.reset_index(drop=True, inplace=True)\n        \n    elif os.path.isfile(data_path_output + name + '.csv'):\n        train_set = pd.read_csv(data_path_output + name + '.csv')\n        train_set.reset_index(drop=True, inplace=True)\n        \n    else:\n        print('training set cannot be found, consider creating it (extract_features = True)')\n\n# Load evaluation set\nname = 'eval_set'\nif name in globals():\n    pass\nelse:\n    if valid == True:\n        name = name + '_valid' \n        \n    if os.path.isfile(data_path_dataset + name + '.csv'):\n        train_set = pd.read_csv(ata_path_dataset + name + '.csv')\n        train_set.reset_index(drop=True, inplace=True)\n        \n    elif os.path.isfile(data_path_output + name + '.csv'):\n        train_set = pd.read_csv(data_path_output + name + '.csv')\n        train_set.reset_index(drop=True, inplace=True)\n        \n    else:\n        print('evaluation set cannot be found, consider creating it (extract_features = True)')","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:50:22.111096Z","iopub.execute_input":"2023-06-10T15:50:22.111813Z","iopub.status.idle":"2023-06-10T15:50:22.123556Z","shell.execute_reply.started":"2023-06-10T15:50:22.111778Z","shell.execute_reply":"2023-06-10T15:50:22.122072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.4.: Train and Store the Individual Regressors ","metadata":{}},{"cell_type":"code","source":"model_description = model_type + ' - Regressor'\nkfold = GroupKFold(N_FOLDS)\ngroup_var = train_set.Subject\ngroups = kfold.split(train_set, groups=group_var)\nregs = []  # List to store trained regressor models\naverage_precision = []  # List to store cross-validation scores\nfprs = []\ntprs = []\nfpr_micros = []\nfpr_macros = []\ntpr_micros = [] \ntpr_macros = [] \nmicro_roc_auc_ovrs = [] \nmacro_roc_auc_ovrs = []\n\n\n# Iterate over the folds\nfor fold, (train_idx, test_idx) in enumerate(tqdm(groups, total=N_FOLDS, desc=\"Folds\")):\n    # Randomly sample a subset of training indices\n    train_idx = pd.Series(train_idx).sample(n=2000000, random_state=42).values  # 2000000\n    #print(tr_idx)\n\n    # Wrap the base regressor with the MultiOutputRegressor\n    multioutput_regressor = CustomMultiOutputRegressor(model)\n    \n\n    # Get training and testing data\n    x_train, y_train = train_set.loc[train_idx, feature_cols].to_numpy(), train_set.loc[train_idx, prediction_cols].to_numpy()\n    x_test, y_test = train_set.loc[test_idx, feature_cols].to_numpy(), train_set.loc[test_idx, prediction_cols].to_numpy()\n\n    # Fit the multioutput regressor\n    multioutput_regressor.fit(x_train, y_train)\n\n    # Append the trained regressor to the list\n    regs.append(multioutput_regressor)\n    pickle.dump(regs, open(model_type + '_regs.sav', 'wb'))\n    \n    y_pred = multioutput_regressor.predict(x_test).clip(0.0, 1.0)\n    \n    fpr, tpr, roc_auc, fpr_micro, fpr_macro, tpr_micro, tpr_macro, micro_roc_auc_ovr, macro_roc_auc_ovr = compute_multiclass_roc_auc_metrics(y_test, y_pred)\n    \n    # Store ROC metrics for each fold\n    fprs.append(fpr)\n    tprs.append(tpr)\n    fpr_micros.append(fpr_micro)\n    fpr_macros.append(fpr_macro)\n    tpr_micros.append(tpr_micro)\n    tpr_macros.append(tpr_macro)\n    micro_roc_auc_ovrs.append(micro_roc_auc_ovr)\n    macro_roc_auc_ovrs.append(macro_roc_auc_ovr)\n    \n    \n    # Evaluate the model on the validation set and calculate average precision score\n    average_precision.append(metrics.average_precision_score(y_test, y_pred))\n\n# Convert the format of the ROC metrics\nfpr_micro_macro = [fpr_micros, fpr_macros]\ndel fpr_micros, fpr_macros\n\ntpr_micro_macro = [tpr_micros, tpr_macros]\ndel tpr_micros, tpr_macros\n\nroc_auc_ovr_micro_macro = [micro_roc_auc_ovrs, macro_roc_auc_ovrs]  \ndel micro_roc_auc_ovrs, macro_roc_auc_ovrs\n     \ngc.collect()\n\n# Evaluate performance\nprint(f'Average precision per fold: {np.round(average_precision,3)}')\nprint(f'Mean average precision: {np.round(np.mean(average_precision), 3)}')\nplot_multiclass_ROC(fprs, tprs, fpr_micro_macro, tpr_micro_macro, roc_auc_ovr_micro_macro)","metadata":{"execution":{"iopub.status.busy":"2023-06-10T15:50:22.125762Z","iopub.execute_input":"2023-06-10T15:50:22.126284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"| Regressor | Mean Average Precision |\n| --- | --- |\n|LGBM | 0.212 |\n|XGB | - |\n|ADAB | - |\n|HGB | 0.201 |\n|RF | - |\n|MLP | - |\n|KN | - |\n|GP| - |\n|XT| - |","metadata":{}},{"cell_type":"markdown","source":"# 3.: Create an Averaging - Ensemble from the Individual Regressors and Evaluate its Performance\n## 3.1.: Create and Store an Ensemble of the Individual Regressors  ","metadata":{}},{"cell_type":"code","source":"# Specify regressors to use for the ensemble\nregressor_ensemble = ['LGBM', 'XGB', 'ADAB',  'HGB'] # 'XT',\n\n# Use the mean average precision of the regressors as averaging weights\nregressor_weights = [0.212, 1, 1, 1]\n\n# Load regressors and create final estimators for the ensemble model\nestimators = []\nfor regressor in regressor_ensemble:\n    match regressor:\n        case 'LGBM':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n        case 'XGB':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n        case 'ADAB':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n        case 'HGB':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n        case 'RF':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n        case 'MLP':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n        case 'KN':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n        case 'GP':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n        case 'XT':\n            estimators.append(pickle.load(open(model_path + regressor +'_reg.sav', 'rb')))\n            \n# Save the ensemble regressors\npickle.dump(estimators, open('ensemble_regressors.sav', 'wb'))\n\n# Save regressor weights\nnp.save('regressor_weights', regressor_weights)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# !!!!!!!!!!!!!!!!!!!!! Include weighted mean (according to mean average precision) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!","metadata":{}},{"cell_type":"markdown","source":"## 3.2.: Evaluate the Performance of the Averaging - Ensemble model","metadata":{}},{"cell_type":"code","source":"# Specify whether to use weighted or unweighted average over the different regressors (weights = mean average precision scores)\nweighted = False\n\nmodel_type = 'Averaging Regressor - Ensemble'\naverage_precision = []  # List to store cross-validation scores\nfprs = []\ntprs = []\nfpr_micros = []\nfpr_macros = []\ntpr_micros = [] \ntpr_macros = [] \nmicro_roc_auc_ovrs = [] \nmacro_roc_auc_ovrs = []\nensemble_predictions = []\naverage_precision = []\n\n# Iterate over the folds\nfor fold in tqdm(range(N_FOLDS)):\n    estimator_predictions = []\n    \n    # Iterate over the different regressors\n    for estimator in estimators:\n        print('!!!!!!!!!! still old ensemble-models instead of new regressor-models dataset !!!!!!!!!!!!!!')\n        estimator_prediction = estimator[fold].predict(eval_set[feature_cols]).clip(0.0, 1.0)\n        estimator_predictions.append(estimator_prediction)\n     \n    # Compute the (weighted) average predictions over the different regressors\n    if weighted == True:\n        ensemble_prediction = np.average(estimator_predictions, axis = 0, weights = regressor_weights)\n        \n    else:\n        ensemble_prediction = np.average(estimator_predictions, axis = 0)\n    \n    fpr, tpr, roc_auc, fpr_micro, fpr_macro, tpr_micro, tpr_macro, micro_roc_auc_ovr, macro_roc_auc_ovr = compute_multiclass_roc_auc_metrics(eval_set[prediction_cols], ensemble_prediction)\n    \n    # Store ROC metrics for each fold\n    fprs.append(fpr)\n    tprs.append(tpr)\n    fpr_micros.append(fpr_micro)\n    fpr_macros.append(fpr_macro)\n    tpr_micros.append(tpr_micro)\n    tpr_macros.append(tpr_macro)\n    micro_roc_auc_ovrs.append(micro_roc_auc_ovr)\n    macro_roc_auc_ovrs.append(macro_roc_auc_ovr)\n    \n    \n    # Evaluate the model on the validation set and calculate average precision score\n    average_precision.append(metrics.average_precision_score(eval_set[prediction_cols], ensemble_prediction))        \n    \n# Convert the format of the ROC metrics\nfpr_micro_macro = [fpr_micros, fpr_macros]\ndel fpr_micros, fpr_macros\n\ntpr_micro_macro = [tpr_micros, tpr_macros]\ndel tpr_micros, tpr_macros\n\nroc_auc_ovr_micro_macro = [micro_roc_auc_ovrs, macro_roc_auc_ovrs]  \ndel micro_roc_auc_ovrs, macro_roc_auc_ovrs\n     \ngc.collect()\n\nprint(f'Average precision per fold: {np.round(average_precision,3)}')\nprint(f'Mean average precision: {np.round(np.mean(average_precision), 3)}')\nplot_multiclass_ROC(fprs, tprs, fpr_micro_macro, tpr_micro_macro, roc_auc_ovr_micro_macro)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}