test_models.py

"""
    Evaluating optimized models with test data
"""

# Libraries
# --------------------------------------------------------------------------------------------------------
import pandas as pd
import numpy as np
from xgboost import XGBClassifier
from sklearn.metrics import confusion_matrix
from sklearn.metrics import f1_score, make_scorer, precision_score, recall_score, accuracy_score, roc_auc_score, average_precision_score
from sklearn.ensemble import RandomForestClassifier, BaggingClassifier, AdaBoostClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.svm import SVC
from sklearn.linear_model import  LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import RocCurveDisplay, roc_curve
from sklearn.metrics import PrecisionRecallDisplay, precision_recall_curve
import matplotlib.pyplot as plt
# --------------------------------------------------------------------------------------------------------

# Reading test data
# --------------------------------------------------------------------------------------------------------
def read_test_data():
    # Load test data
    X_test_pre = np.load('../gen_train_data/data/output/pre/X_test_pre.npy', allow_pickle=True)
    y_test_pre = np.load('../gen_train_data/data/output/pre/y_test_pre.npy', allow_pickle=True)
    X_test_post = np.load('../gen_train_data/data/output/post/X_test_post.npy', allow_pickle=True)
    y_test_post = np.load('../gen_train_data/data/output/post/y_test_post.npy', allow_pickle=True)

    # Load ORIGINAL training data
    X_train_pre = np.load('../gen_train_data/data/output/pre/X_train_pre.npy', allow_pickle=True)
    y_train_pre = np.load('../gen_train_data/data/output/pre/y_train_pre.npy', allow_pickle=True)
    X_train_post = np.load('../gen_train_data/data/output/post/X_train_post.npy', allow_pickle=True)
    y_train_post = np.load('../gen_train_data/data/output/post/y_train_post.npy', allow_pickle=True)

    # Load oversampled training data
    X_train_over_pre = np.load('../gen_train_data/data/output/pre/X_train_over_pre.npy', allow_pickle=True)
    y_train_over_pre = np.load('../gen_train_data/data/output/pre/y_train_over_pre.npy', allow_pickle=True)
    X_train_over_post = np.load('../gen_train_data/data/output/post/X_train_over_post.npy', allow_pickle=True)
    y_train_over_post = np.load('../gen_train_data/data/output/post/y_train_over_post.npy', allow_pickle=True)

    # Load undersampled training data
    X_train_under_pre = np.load('../gen_train_data/data/output/pre/X_train_under_pre.npy', allow_pickle=True)
    y_train_under_pre = np.load('../gen_train_data/data/output/pre/y_train_under_pre.npy', allow_pickle=True)
    X_train_under_post = np.load('../gen_train_data/data/output/post/X_train_under_post.npy', allow_pickle=True)
    y_train_under_post = np.load('../gen_train_data/data/output/post/y_train_under_post.npy', allow_pickle=True)

    data_dic = {
        "X_test_pre": X_test_pre,
        "y_test_pre": y_test_pre,
        "X_test_post": X_test_post,
        "y_test_post": y_test_post,
        "X_train_pre": X_train_pre,
        "y_train_pre": y_train_pre,
        "X_train_post": X_train_post,
        "y_train_post": y_train_post,
        "X_train_over_pre": X_train_over_pre,
        "y_train_over_pre": y_train_over_pre,
        "X_train_over_post": X_train_over_post,
        "y_train_over_post": y_train_over_post,
        "X_train_under_pre": X_train_under_pre,
        "y_train_under_pre": y_train_under_pre,
        "X_train_under_post": X_train_under_post,
        "y_train_under_post": y_train_under_post,
    }

    return data_dic
# --------------------------------------------------------------------------------------------------------

# Returning tuned models for each situation
# --------------------------------------------------------------------------------------------------------
def get_tuned_models(group_id, method_id):
    # 1. PRE
    if group_id == 0:
        # 1.1) Trained with original dataset
        if method_id == 0:
            tuned_models = {
            "DT" : DecisionTreeClassifier(**{'splitter': 'best', 'max_features': 'sqrt', 'criterion': 'entropy'}), 
            "RF" : RandomForestClassifier(**{'criterion': 'entropy', 'max_features': 'sqrt', 'n_estimators': 123}), 
            "Bagging" : BaggingClassifier(**{'max_features': 1.0, 'max_samples': 0.8, 'n_estimators': 13, 'warm_start': False}),
            "AB" : AdaBoostClassifier(**{'learning_rate': 1.8473150336970519, 'n_estimators': 96, 'algorithm': 'SAMME'}), 
            "XGB": XGBClassifier(**{'learning_rate': 0.21528982071549305, 'max_depth': 6, 'n_estimators': 804}),
            "LR" : LogisticRegression(**{'solver': 'lbfgs', 'penalty': 'l2','max_iter': 1000}), 
            "SVM" : SVC(**{'C': 1.051871311397777, 'kernel': 'linear', 'max_iter':1000, 'probability': True}), 
            "MLP" : MLPClassifier(**{'activation': 'identity', 'hidden_layer_sizes': 78, 'learning_rate': 'constant','max_iter':500})
            }
        # 1.2) Trained with original dataset and cost-sensitive learning
        elif method_id == 1:
            tuned_models = {
            "DT": DecisionTreeClassifier(**{'splitter': 'best', 'max_features': 'log2', 'criterion': 'entropy', 'class_weight': 'balanced'}),
            "RF": RandomForestClassifier(**{'criterion': 'entropy', 'max_features': 'sqrt', 'n_estimators': 238, 'class_weight': 'balanced'}),
            "Bagging": BaggingClassifier(**{'max_features': 1.0, 'max_samples': 0.8, 'n_estimators': 22, 'warm_start': False, 'estimator': DecisionTreeClassifier(class_weight='balanced')}),
            "AB": AdaBoostClassifier(**{'learning_rate': 1.7136783954287846, 'n_estimators': 99, 'algorithm': 'SAMME', 'estimator': DecisionTreeClassifier(class_weight='balanced')}),
            "LR": LogisticRegression(**{'solver': 'lbfgs', 'penalty': 'l2', 'max_iter': 1000, 'class_weight': 'balanced'}),
            "SVM": SVC(**{'C': 1.480857958217729, 'kernel': 'linear', 'max_iter': 1000, 'class_weight': 'balanced', 'probability': True}),
            }
        # 1.3) Trained with oversampled training dataset
        elif method_id == 2:
            tuned_models = {
            "DT" : DecisionTreeClassifier(**{'splitter': 'best', 'max_features': 'sqrt', 'criterion': 'log_loss'}), 
            "RF" : RandomForestClassifier(**{'criterion': 'gini', 'max_features': 'sqrt', 'n_estimators': 121}), 
            "Bagging" : BaggingClassifier(**{'max_features': 1.0, 'max_samples': 1.0, 'n_estimators': 22, 'warm_start': True}),
            "AB" : AdaBoostClassifier(**{'learning_rate': 1.4640913091426446, 'n_estimators': 145, 'algorithm': 'SAMME'}), 
            "XGB": XGBClassifier(**{'learning_rate': 0.19621698151985992, 'max_depth': 7, 'n_estimators': 840}),
            "LR" : LogisticRegression(**{'solver': 'lbfgs', 'penalty': 'l2', 'max_iter': 1000}), 
            "SVM" : SVC(**{'C': 1.590799972846728, 'kernel': 'poly', 'max_iter':1000, 'probability': True}), 
            "MLP" : MLPClassifier(**{'activation': 'relu', 'hidden_layer_sizes': 112, 'learning_rate': 'constant', 'max_iter':500})
            }
        # 1.4) Trained with undersampled training dataset
        elif method_id == 3:
            tuned_models = {
            "DT" : DecisionTreeClassifier(**{'splitter': 'best', 'max_features': 'sqrt', 'criterion': 'log_loss'}), 
            "RF" : RandomForestClassifier(**{'criterion': 'gini', 'max_features': 'sqrt', 'n_estimators': 148}), 
            "Bagging" : BaggingClassifier(**{'max_features': 1.0, 'max_samples': 0.8, 'n_estimators': 24, 'warm_start': True}),
            "AB" : AdaBoostClassifier(**{'learning_rate': 1.7970533619575801, 'n_estimators': 122, 'algorithm': 'SAMME'}), 
            "XGB": XGBClassifier(**{'learning_rate': 0.13148624656904934, 'max_depth': 9, 'n_estimators': 723}),
            "LR" : LogisticRegression(**{'solver': 'sag', 'penalty': 'l2', 'max_iter': 1000}), 
            "SVM" : SVC(**{'C': 1.383651513577477, 'kernel': 'poly', 'max_iter':1000, 'probability': True}), 
            "MLP" : MLPClassifier(**{'activation': 'relu', 'hidden_layer_sizes': 89, 'learning_rate': 'invscaling', 'max_iter':500})
            }
    # 2. POST
    else:
        # 2.1) Trained with original dataset
        if method_id == 0:
            tuned_models = {
            "DT" : DecisionTreeClassifier(**{'splitter': 'best', 'max_features': 'sqrt', 'criterion': 'log_loss'}), 
            "RF" : RandomForestClassifier(**{'criterion': 'entropy', 'max_features': 'sqrt', 'n_estimators': 120}), 
            "Bagging" : BaggingClassifier(**{'max_features': 1.0, 'max_samples': 0.8, 'n_estimators': 38, 'warm_start': True}),
            "AB" : AdaBoostClassifier(**{'learning_rate': 1.9069394544838472, 'n_estimators': 121, 'algorithm': 'SAMME'}), 
            "XGB": XGBClassifier(**{'learning_rate': 0.24787889985627387, 'max_depth': 4, 'n_estimators': 956}),
            "LR" : LogisticRegression(**{'solver': 'lbfgs', 'penalty': 'l2'}), 
            "SVM" : SVC(**{'C': 1.7965537393241109, 'kernel': 'linear', 'max_iter':1000, 'probability': True}), 
            "MLP" : MLPClassifier(**{'activation': 'relu', 'hidden_layer_sizes': 147, 'learning_rate': 'invscaling', 'max_iter':500})
            }
        # 2.2) Trained with original dataset and cost-sensitive learning
        elif method_id == 1:
            tuned_models = {
            "DT": DecisionTreeClassifier(**{'splitter': 'best', 'max_features': 'sqrt', 'criterion': 'gini', 'class_weight': 'balanced'}),
            "RF": RandomForestClassifier(**{'criterion': 'entropy', 'max_features': 'sqrt', 'n_estimators': 138, 'class_weight': 'balanced'}),
            "Bagging": BaggingClassifier(**{'max_features': 1.0, 'max_samples': 1.0, 'n_estimators': 66, 'warm_start': True, 'estimator': DecisionTreeClassifier(class_weight='balanced')}),
            "AB": AdaBoostClassifier(**{'learning_rate': 1.92541653518023, 'n_estimators': 114, 'algorithm': 'SAMME', 'estimator': DecisionTreeClassifier(class_weight='balanced')}),
            "LR": LogisticRegression(**{'solver': 'lbfgs', 'penalty': 'l2', 'max_iter': 1000, 'class_weight': 'balanced'}),
            "SVM": SVC(**{'C': 0.8395104850983046, 'kernel': 'linear', 'max_iter': 1000, 'class_weight': 'balanced', 'probability': True})
            }
        # 2.3) Trained with oversampled training dataset
        elif method_id == 2:
            tuned_models = {
            "DT" : DecisionTreeClassifier(**{'splitter': 'best', 'max_features': 'log2', 'criterion': 'entropy'}), 
            "RF" : RandomForestClassifier(**{'criterion': 'gini', 'max_features': 'sqrt', 'n_estimators': 118}), 
            "Bagging" : BaggingClassifier(**{'max_features': 1.0, 'max_samples': 1.0, 'n_estimators': 56, 'warm_start': False}),
            "AB" : AdaBoostClassifier(**{'learning_rate': 1.5933610622176648, 'n_estimators': 114, 'algorithm': 'SAMME'}), 
            "XGB": XGBClassifier(**{'learning_rate': 0.059934879882855396, 'max_depth': 9, 'n_estimators': 660}),
            "LR" : LogisticRegression(**{'solver': 'lbfgs', 'penalty': 'l2', 'max_iter': 1000}), 
            "SVM" : SVC(**{'C': 1.2237930722499044, 'kernel': 'poly', 'max_iter':1000, 'probability': True}), 
            "MLP" : MLPClassifier(**{'activation': 'identity', 'hidden_layer_sizes': 134, 'learning_rate': 'invscaling', 'max_iter':500})
            }
        # 2.4) Trained with undersampled training dataset
        elif method_id == 3:
            tuned_models = {
            "DT" : DecisionTreeClassifier(**{'splitter': 'best', 'max_features': 'log2', 'criterion': 'log_loss'}), 
            "RF" : RandomForestClassifier(**{'criterion': 'gini', 'max_features': 'sqrt', 'n_estimators': 151}), 
            "Bagging" : BaggingClassifier(**{'max_features': 1.0, 'max_samples': 1.0, 'n_estimators': 20, 'warm_start': False}),
            "AB" : AdaBoostClassifier(**{'learning_rate': 1.6523810056317618, 'n_estimators': 89, 'algorithm': 'SAMME'}), 
            "XGB": XGBClassifier(**{'learning_rate': 0.18430397856234193, 'max_depth': 4, 'n_estimators': 956}),
            "LR" : LogisticRegression(**{'solver': 'lbfgs', 'penalty': 'l2', 'max_iter': 1000}), 
            "SVM" : SVC(**{'C': 1.1807459108651588, 'kernel': 'linear', 'max_iter':1000, 'probability': True}), 
            "MLP" : MLPClassifier(**{'activation': 'identity', 'hidden_layer_sizes': 55, 'learning_rate': 'constant', 'max_iter':500})
            }
    return tuned_models
# --------------------------------------------------------------------------------------------------------

# Scorers
# --------------------------------------------------------------------------------------------------------
def TN_scorer(clf, X, y):
    """Gives the number of samples predicted as true negatives"""
    y_pred = clf.predict(X)
    cm = confusion_matrix(y, y_pred)
    TN = cm[0,0]
    return TN
def FN_scorer(clf, X, y):
    """Gives the number of samples predicted as false negatives"""
    y_pred = clf.predict(X)
    cm = confusion_matrix(y, y_pred)
    FN = cm[0,1]
    return FN
def FP_scorer(clf, X, y):
    """Gives the number of samples predicted as false positive"""
    y_pred = clf.predict(X)
    cm = confusion_matrix(y, y_pred)
    FP = cm[1,0]
    return FP
def TP_scorer(clf, X, y):
    """Gives the number of samples predicted as true positive"""
    y_pred = clf.predict(X)
    cm = confusion_matrix(y, y_pred)
    TP = cm[1,1]
    return TP

def negative_recall_scorer(clf, X, y):
    """Gives the negative recall defined as the (number of true_negative_samples)/(total number of negative samples)"""
    y_pred = clf.predict(X)
    cm = confusion_matrix(y, y_pred)
    TN_prop = cm[0,0]/(cm[0,1]+cm[0,0])
    return TN_prop
# --------------------------------------------------------------------------------------------------------

if __name__ == "__main__":
    # Reading testing data
    data_dic = read_test_data()

    # Setup
    # --------------------------------------------------------------------------------------------------------
    # Scorings to use for model evaluation
    scorings = {
        'F1':make_scorer(f1_score), 
        'NREC': negative_recall_scorer, 
        'REC':make_scorer(recall_score), 
        'PREC':make_scorer(precision_score), 
        'ACC': make_scorer(accuracy_score),
        'TN':TN_scorer, 
        'FN':FN_scorer, 
        'FP':FP_scorer, 
        'TP':TP_scorer,
        'AUROC': make_scorer(roc_auc_score),  # AUROC requires decision function or probability outputs
        'AUPRC': make_scorer(average_precision_score)  # AUPRC requires probability outputs
        } 
    method_names = {
        0: "ORIG",
        1: "ORIG_CW",
        2: "OVER",
        3: "UNDER"
    }
    # --------------------------------------------------------------------------------------------------------

    # Evaluating performance using test dataset
    # --------------------------------------------------------------------------------------------------------
    scores_sheets = {} # To store score dfs as sheets in the same excel file
    for i, group in enumerate(['pre', 'post']):
        # Get test dataset based on group
        X_test = data_dic['X_test_' + group]
        y_test = data_dic['y_test_' + group]
        for j, method in enumerate(['', '', 'over_', 'under_']):
            print(f"{group}-{method}")
            # Get train dataset based on group and method
            X_train = data_dic['X_train_' + method + group]
            y_train = data_dic['y_train_' + method + group]
            # Get tuned models for this group and method
            models = get_tuned_models(group_id=i, method_id=j)
            # Scores df
            scores_df = pd.DataFrame(index=models.keys(), columns=scorings.keys())
            # Create a figure for all models in this group-method
            fig, axes = plt.subplots(len(models), 2, figsize=(12, 8 * len(models)))
            if len(models) == 1:  # Adjustment if there's only one model (axes indexing issue)
                axes = [axes]
            # Evaluate each model
            for model_idx, (model_name, model) in enumerate(models.items()):
                # ----------- TEMPORAL -------------
                if model_name == "DT" or model_name == "RF":
                    # Train the model (it was just initialized above)
                    model.fit(X_train, y_train)
                    if hasattr(model, "decision_function"):
                        y_score = model.decision_function(X_test)
                    else:
                        y_score = model.predict_proba(X_test)[:, 1]  # Use probability of positive class
                    # Calculate ROC curve and ROC area for each class
                    fpr, tpr, _ = roc_curve(y_test, y_score, pos_label=model.classes_[1])
                    roc_display = RocCurveDisplay(fpr=fpr, tpr=tpr).plot(ax=axes[model_idx][0]) 
                    # Calculate precision-recall curve
                    precision, recall, _ = precision_recall_curve(y_test, y_score, pos_label=model.classes_[1])
                    pr_display = PrecisionRecallDisplay(precision=precision, recall=recall).plot(ax=axes[model_idx][1])
                    axes[model_idx][0].set_title(f'ROC Curve for {model_name}')
                    axes[model_idx][1].set_title(f'PR Curve for {model_name}')
                    # Evaluate at each of the scores of interest
                    for score_name, scorer in scorings.items():
                        score_value = scorer(model, X_test, y_test)
                        scores_df.at[model_name, score_name] = score_value
            # Adjust layout and save/show figure
            plt.tight_layout()
            plt.savefig(f'./test_results/roc_pr_curves/{group}_{method_names[j]}.svg', format='svg', dpi=500)
            plt.close(fig)
            # Store the DataFrame in the dictionary with a unique key for each sheet
            sheet_name = f"{group}_{method_names[j]}"
            scores_sheets[sheet_name] = scores_df
    # Write results to Excel file
    with pd.ExcelWriter('./test_results/testing_tuned_models.xlsx') as writer:
        for sheet_name, data in scores_sheets.items():
            data.to_excel(writer, sheet_name=sheet_name)
    # --------------------------------------------------------------------------------------------------------