classifier.py

d#!interpreter [optional-arg]
# -*- coding: utf-8 -*-

__author__ = ["David Cardona-Vasquez"]
__copyright__ = "Copyright 2023-2025, Graz University of Technology"
__credits__ = ["David Cardona-Vasquez"]
__license__ = "MIT"
__maintainer__ = "David Cardona-Vasquez"
__status__ = "Development"

import matplotlib.pyplot as plt
import seaborn as sns

# the data module is used to load the data from the csv files and data wrangling for the models
from data import *


from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, confusion_matrix, classification_report, balanced_accuracy_score
from sklearn.model_selection import train_test_split


def create_input_classifier(case:str, sep=';', sep_dec=',', header=None) -> pd.DataFrame:
    """This function creates all the candidate features for the classifier.
    The function loads the data from the case's csv files and performs the necessary data transformations
     to create the features."""

    renewables_data = load_renewables(case, sep_dec, sep)
    cf_data = load_cf(case, sep_dec, sep, header=0)
    demand_data = load_demand(case, sep_dec, sep, header=0)
    network_data = load_network(case, sep_dec)

    df_input = cf_data[0]['cf'].copy()
    df_input['resource'] = cf_data[0]['generator']
    df_input['max_gen'] = df_input['cf'] * \
                          renewables_data.loc[renewables_data.unit == cf_data[0]['generator']]['max'].values[0]
    for x in cf_data[1:]:
        y = x['cf'].copy()
        y['resource'] = x['generator']
        y['max_gen'] = y['cf'] * renewables_data.loc[renewables_data.unit == x['generator']]['max'].values[0]
        df_input = pd.concat([df_input, y])

    df_input.reset_index(drop=True, inplace=True)

    df_demand = demand_data[0]['demand'].copy()
    df_demand['bus'] = demand_data[0]['bus']
    for x in demand_data[1:]:
        y = x['demand'].copy()
        y['bus'] = x['bus']
        df_demand = pd.concat([df_demand, y])

    df_input = pd.merge(df_input, renewables_data, left_on='resource', right_on='unit', how='left')
    df_input.drop(columns=['max', 'min', 'type', 'unit'], inplace=True)

    df_input_agg = df_input.groupby(by=['hour', 'bus']).agg({'max_gen': 'sum'}).reset_index()

    df_exp_cap = network_data.groupby(by='from').sum().reset_index()
    if 'line' in network_data.columns:
        df_exp_cap.drop(columns=['line'], inplace=True)
    df_exp_cap.drop(columns=['to'], inplace=True)
    df_exp_cap.columns = ['bus', 'exp_cap']
    df_input_agg = pd.merge(df_input_agg, df_exp_cap, left_on='bus', right_on='bus', how='left')
    df_input_agg['max_gen_exp_cap_ratio'] = df_input_agg.max_gen / df_input_agg.exp_cap
    df_input_agg['max_gen_exp_cap_ratio_lag_1'] = df_input_agg.max_gen.shift(1) / df_input_agg.exp_cap.shift(1)
    df_input_agg['max_gen_exp_cap_ratio_lead_1'] = df_input_agg.max_gen.shift(-1) / df_input_agg.exp_cap.shift(-1)
    df_input_agg.drop(columns=['exp_cap'], inplace=True)

    df_demand.reset_index(drop=True, inplace=True)
    df_demand_agg = df_demand.groupby(by=['hour']).sum('demand').reset_index()

    df_input_agg = df_input_agg.groupby(by='hour').agg({'max_gen_exp_cap_ratio': 'mean',
                                                        'max_gen': 'sum', 'max_gen_exp_cap_ratio_lag_1': 'mean',
                                                        'max_gen_exp_cap_ratio_lead_1': 'mean'}).reset_index()

    ###### Start of feature creation process
    #---------------------------------------
    df_classifier_input = pd.merge(df_input_agg, df_demand_agg, on=['hour'], how='left')
    df_classifier_input['net_demand'] = df_classifier_input.demand - df_classifier_input.max_gen
    df_classifier_input['net_demand_lag_1'] = df_classifier_input.net_demand.shift(1)
    df_classifier_input['net_demand_lag_2'] = df_classifier_input.net_demand_lag_1.shift(1)
    df_classifier_input['net_demand_lag_3'] = df_classifier_input.net_demand_lag_2.shift(1)
    df_classifier_input['net_demand_lead_1'] = df_classifier_input.net_demand.shift(-1)
    df_classifier_input['net_demand_lead_2'] = df_classifier_input.net_demand_lead_1.shift(-1)
    df_classifier_input['net_demand_lead_3'] = df_classifier_input.net_demand_lead_2.shift(-1)

    df_classifier_input['demand_lag_1'] = df_classifier_input.demand.shift(1)
    df_classifier_input['demand_lag_2'] = df_classifier_input.demand_lag_1.shift(1)
    df_classifier_input['demand_lag_3'] = df_classifier_input.demand_lag_2.shift(1)
    df_classifier_input['demand_lead_1'] = df_classifier_input.demand.shift(-1)
    df_classifier_input['demand_lead_2'] = df_classifier_input.demand_lead_1.shift(-1)
    df_classifier_input['demand_lead_3'] = df_classifier_input.demand_lead_2.shift(-1)
    df_classifier_input['max_gen_lag_1'] = df_classifier_input.max_gen.shift(1)
    df_classifier_input['max_gen_lag_2'] = df_classifier_input.max_gen_lag_1.shift(1)
    df_classifier_input['max_gen_lag_3'] = df_classifier_input.max_gen_lag_2.shift(1)
    df_classifier_input['max_gen_lead_1'] = df_classifier_input.max_gen.shift(-1)
    df_classifier_input['max_gen_lead_2'] = df_classifier_input.max_gen_lead_1.shift(-1)
    df_classifier_input['max_gen_lead_3'] = df_classifier_input.max_gen_lead_2.shift(-1)

    df_classifier_input['ren_demand_ratio'] = df_classifier_input.max_gen / df_classifier_input.demand
    df_classifier_input['ren_demand_ratio_lag_1'] = df_classifier_input.max_gen_lag_1 / df_classifier_input.demand_lag_1
    df_classifier_input['ren_demand_ratio_lag_2'] = df_classifier_input.max_gen_lag_2 / df_classifier_input.demand_lag_2
    df_classifier_input['ren_demand_ratio_lag_3'] = df_classifier_input.max_gen_lag_3 / df_classifier_input.demand_lag_3
    df_classifier_input['ren_demand_ratio_lead_1'] = df_classifier_input.max_gen_lead_1 / df_classifier_input.demand_lead_1
    df_classifier_input['ren_demand_ratio_lead_2'] = df_classifier_input.max_gen_lead_2 / df_classifier_input.demand_lead_2
    df_classifier_input['ren_demand_ratio_lead_3'] = df_classifier_input.max_gen_lead_3 / df_classifier_input.demand_lead_3

    ######### End of feature creation process
    # ---------------------------------------

    df_classifier_input.set_index('hour', inplace=True)

    return df_classifier_input


def determine_feature_importance(X: pd.DataFrame, y:pd.Series) -> np.ndarray:
    """
    This function determines the feature importance of the features in X using a logistic regression model.
    :param X: input features
    :param y: target data
    :return: importances based on the coefficients of the logistic regression model
    """

    scaler = MinMaxScaler()  # Scale the data for logistic regression
    X_scaled = scaler.fit_transform(X)

    # Logistic Regression
    logreg = LogisticRegression(max_iter=10000)
    logreg.fit(X_scaled, y)
    coefficients = logreg.coef_[0]
    feature_importance_logreg = pd.Series(coefficients, index=X.columns)
    feature_importance_logreg = feature_importance_logreg.abs()

    feature_importance_logreg = feature_importance_logreg / feature_importance_logreg.sum()

    return feature_importance_logreg

def fit_model(X_train: pd.DataFrame, y_train:pd.Series) -> RandomForestClassifier:
    """
    Fits a random forest classifier to the training data and returns the fitted model.
    :param X_train: input features
    :param y_train: target data
    :return: a fitted random forest classifier
    """
    rf = RandomForestClassifier(random_state=42, n_estimators=10000, n_jobs=-1)
    rf.fit(X_train, y_train)

    return rf

def prepare_classifier_input(df_classifier_input: pd.DataFrame, features: np.ndarray, threshold=0.025, test_size=0.2) -> tuple:
    """
    Performs the train-test split of the data and returns the training and test sets.
    :param df_classifier_input: candidate features to fit the classifier
    :param features: vector of feature importance in percentage
    :param threshold: cut-off value for the feature importance to consider in the classifier
    :param test_size: proportion of the dataset to include in the test split
    :return: training and test sets for the classifier with their respective target values
    """

    df_input =df_classifier_input.iloc[:, :-1]
    df_target = df_classifier_input.iloc[:, -1]
    idx_features = features > threshold
    X = df_input.loc[:, idx_features].values
    y = df_target.values


    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, stratify=y, random_state=42)


    return X_train, X_test, y_train, y_test


if __name__ == '__main__':

    # feature importance to consider in the classifier
    threshold = 0.04

    # 2023 is the test year
    df_classifier_input_2023 = create_input_classifier('austria_2023', sep=';', sep_dec='.')
    target_2023 = pd.read_csv('data/austria_2023/decomposition.csv', sep=';')
    df_classifier_input_2023 = df_classifier_input_2023.merge(target_2023, left_index=True, right_on='period', how='left')
    df_classifier_input_2023.set_index('period', inplace=True)

    # 2022 is the training year
    df_classifier_input_2022 = create_input_classifier('austria_2022', sep=';', sep_dec='.')
    target_2022 = pd.read_csv('data/austria_2022/decomposition.csv', sep=';')
    df_classifier_input_2022 = df_classifier_input_2022.merge(target_2022, left_index=True, right_on='period', how='left')
    df_classifier_input_2022.set_index('period', inplace=True)


    df_classifier_input_2022.dropna(inplace=True)
    features = determine_feature_importance(df_classifier_input_2022.iloc[:, :-1], df_classifier_input_2022.iloc[:, -1])
    X_train, X_test, y_train, y_test = prepare_classifier_input(df_classifier_input_2022, features,
                                                                threshold=threshold, test_size=0.2)
    rf = fit_model(X_train, y_train)
    y_pred_rf_new = rf.predict(X_test)

    # replicating some of the plots from the manuscript
    from matplotlib.ticker import PercentFormatter

    plt.figure(figsize=(62, 35))
    features = features.sort_values(ascending=False)
    sns.barplot(x=features.values, y=features.index, color='darkblue')
    plt.axvline(x=threshold, color='r', linestyle='--', linewidth=3)
    plt.title('Feature importance from logistic regression', fontsize=45,
              fontweight='bold', pad=30)
    plt.yticks(fontsize=35, fontstyle='italic')
    plt.xticks(fontsize=30)
    plt.gca().xaxis.set_major_formatter(PercentFormatter(1))
    plt.grid(True)
    plt.show()



    print("Accuracy:", accuracy_score(y_test, y_pred_rf_new))
    print("Accuracy (Balanced):", balanced_accuracy_score(y_test, y_pred_rf_new))
    print("\nConfusion Matrix:\n", confusion_matrix(y_test, y_pred_rf_new))
    print("\nClassification Report:\n", classification_report(y_test, y_pred_rf_new))

    df_aux = df_classifier_input_2023.iloc[:, :-1]
    df_aux = df_aux.loc[:, features > threshold]
    y_pred_2023 = rf.predict(df_aux.values)
    y_true_2023 = df_classifier_input_2023.iloc[:, -1].values

    print("Accuracy:", accuracy_score(y_true_2023, y_pred_2023))
    print("Accuracy (Balanced):", balanced_accuracy_score(y_true_2023, y_pred_2023))
    print("\nConfusion Matrix:\n", confusion_matrix(y_true_2023, y_pred_2023))
    print("\nClassification Report:\n", classification_report(y_true_2023, y_pred_2023))

    y_pred_2023.drop(columns=[y_pred_2023.columns[0]],inplace=True)