<reponame>Songyosk/GBSFS4MPP
"""
Module to train, optimize and evulate ML model 

Author: <NAME>
Email: <EMAIL>
"""  

import os   
import numpy as np                                                                                                                                                                               
import pandas as pd
import joblib
import matplotlib.pyplot as plt

import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from tensorflow.keras.layers import BatchNormalization, Dropout
from tensorflow.python.keras import backend as K
from tensorflow.keras.optimizers import Adam, SGD, Adadelta, RMSprop
from tensorflow.keras.optimizers.schedules import ExponentialDecay
from tensorflow.keras.callbacks import EarlyStopping
from tensorflow.keras import metrics
from tensorflow.keras.callbacks import ModelCheckpoint
from tensorflow.keras.utils import to_categorical

from sklearn.preprocessing import label_binarize
from sklearn.metrics import multilabel_confusion_matrix, roc_curve, roc_auc_score, max_error, \
                            auc, f1_score, classification_report, recall_score, precision_recall_curve, \
                            balanced_accuracy_score, confusion_matrix, accuracy_score, average_precision_score, \
                            hamming_loss, matthews_corrcoef, mean_squared_error, mean_absolute_error, r2_score, \
                            plot_confusion_matrix, explained_variance_score


from lightgbm.sklearn import LGBMClassifier, LGBMRegressor
from xgboost import XGBClassifier, XGBRegressor

from skopt import forest_minimize, gbrt_minimize, gp_minimize, dummy_minimize
from sklearn.model_selection import cross_val_score
from skopt.utils import use_named_args
from skopt.space import Real, Integer
from skopt import dump, load
from skopt.plots import plot_convergence, plot_objective, plot_evaluations
from skopt import dump, load

import statsmodels.api as sm
import statsmodels.formula.api as smf

from itertools import cycle


class optimization():
    """
    Optimize and evaluate ML model 

    args: 
        (1) path_to_train_data (type:str) - location of the training data 
        (2) path_to_test_data (type:str) - location of the test data 
        (3) path_to_features (type:str) - location of the features to use
        (4) path_to_save (type:str) - location to save new data files
        (5) problem (type:str) - whether it is a 'classification' or 'regression' problem
    
    return: 
        (1) performance evaluation of ML model
    """

    def __init__(self, path_to_train_data, path_to_test_data, path_to_features, path_to_save, problem, *args, **kwargs):
        self.path_to_save = path_to_save
        self.sample_train = joblib.load(path_to_train_data)  
        self.sample_test = joblib.load(path_to_test_data)

        self.sample_train.columns = self.sample_test.columns
        self.RFE_features = joblib.load(path_to_features)

        self.features = self.sample_train.columns.values[:-1]
        self.target = self.sample_train.columns.values[-1]

        print('Name of target column: ', self.target)
        print('No. of exploratory features: ', len(self.RFE_features))

        self.problem = problem

        self.target_classes = kwargs.get('target_classes')
        self.multiclass = kwargs.get('multiclass')



    def set_hyperparameters(self):
        """
        Define the hyperparameter space where optimization will be conducted
        """

        dim_learning_rate = Real(low = 1e-7, high = 1e-3, prior = 'log-uniform', name = 'learning_rate')
        dim_num_dense_layers = Integer(low = 1, high = 4, name = 'num_dense_layers')
        dim_num_dense_nodes = Integer(low = 50, high = 500, name = 'num_dense_nodes')
        dim_batch_size = Integer(low = 1, high = 4, name = 'batch_size')


        self.dimensions = [
                            dim_learning_rate,
                            dim_num_dense_layers,
                            dim_num_dense_nodes,
                            #dim_activation,
                            dim_batch_size
                            ]

        self.hyperparameters = [
                                'learning_rate',
                                'num_dense_layers',
                                'num_dense_nodes',
                                'batch_size'
                                ]

        # Initial guesses
        lr_0 = 1e-04
        dl_0 = 2
        dn_0 = 200
        batch = 2

        self.default_parameters = [lr_0, dl_0, dn_0, batch]


        return self.dimensions, self.default_parameters



    def set_callbacks(self):
        """
        Define early stoppting criteria to be used during the training / optimistation process
        """

        if self.problem == 'regression':
            self.callback = [
                            EarlyStopping(
                                        monitor = 'val_loss', 
                                        patience = 5,
                                        min_delta = 0,
                                        verbose = 1, 
                                        mode = 'min',
                                        restore_best_weights = True
                                        )]


        elif self.problem == 'classification':
            self.callback = [
                            EarlyStopping(
                                        monitor = 'val_loss', 
                                        patience = 5,
                                        min_delta = 0,
                                        verbose = 1, 
                                        mode = 'max',
                                        restore_best_weights = True
                                        )]


        return self.callback



    def create_model(self, learning_rate, num_dense_layers, num_dense_nodes, batch_size):
        """
        Create ANN architecture for either a classificaitonn or a regression problem, with the specified performance metrics.
        This function is called under the 'run()' function

        args:
            (1) learning_rate (type:float) - the learning rate for the model
            (2) num_dense_layers (type:int) - the number of dense layers in the model
            (3) num_dense_nodes (type:int) - the number of dense nodes in the layers
            (4) batch_size (type:int) - the batch size used to train the model (the number is multiplied by 8)

        return: 
            (1) ANN model
        """

        # Dimension of first layer
        input_dim = self.sample_train[self.RFE_features].shape[1]


        # Define activation to use 
        if self.problem == 'classification':
            activation = 'sigmoid' 

        elif self.problem == 'regression':
            activation = 'relu'


        # Start the model making process and create our first layer
        model = Sequential()
        model.add(Dense(units = input_dim, 
                        input_dim = input_dim,
                        kernel_initializer = 'normal',
                        #bias_initializer='zeros',
                        activation = activation
                       ))
        

        # Create a loop making a new dense layer for the amount passed to this model
        for i in range(num_dense_layers):
            name = 'layer_dense_{0}'.format(i+1)
            model.add(Dense(num_dense_nodes, 
                            kernel_initializer = 'normal',
                            #bias_initializer = 'zeros',
                            activation = activation
                            ))
        
        # Create last layer
        if self.problem == 'regression':
            # Add regression layer.
            model.add(Dense(units = 1, kernel_initializer = 'normal', activation = 'relu'))


        elif self.problem == 'classification':
            if self.multiclass == True:
                model.add(Dense(units = self.target_classes, activation='softmax'))
            
            elif self.multiclass == False:
                mode.add(Dense(units = 1, activation = 'softmax'))


        # Setup optimizer 
        lr_schedule = ExponentialDecay(initial_learning_rate = learning_rate,
                                       decay_steps = 100000,
                                       decay_rate = 0.96
                                       )

        #optimizer = SGD(learning_rate=lr_schedule)
        #optimizer = RMSprop(learning_rate=lr_schedule, momentum=0.9)
        #optimizer = Adadelta(learning_rate=lr_schedule)
        optimizer = Adam(learning_rate = lr_schedule)


        # Compile using optimizer 
        if self.problem == 'regression':
            model.compile(
                        optimizer = optimizer,
                        loss = 'mean_squared_error',
                        metrics = [
                                    metrics.MeanSquaredError(name = 'mean_squared_error'),
                                    metrics.RootMeanSquaredError(name = 'root_mean_squared_error'),
                                    metrics.MeanAbsoluteError(name = 'mean_absolute_error')
                                    ])


        elif self.problem == 'classification':
            if self.multiclass == True:
                model.compile(
                            loss = 'categorical_crossentropy', 
                            optimizer = optimizer, 
                            metrics = [
                                        metrics.Accuracy(name = 'accuracy'),
                                        metrics.Precision(name = 'precision'),
                                        metrics.Recall(name = 'recall')
                                        ])


            elif self.multiclass == False:
            model.compile(
                        loss = 'binary_crossentropy', 
                        optimizer = optimizer, 
                        metrics = [
                                    metrics.Accuracy(name = 'accuracy'),
                                    metrics.Precision(name = 'precision'),
                                    metrics.Recall(name = 'recall')
                                    ])


        return model



    def run(self, optimization_method):
        """
        Execute optimization using one of the methods:
        (1) 'random_search' -> dummy_minimize
        (2) 'bayesian' -> gp_minimize, 
        (3) 'gradient_bossted_trees' -> gbrt_minimize, 
        (4) 'decision_trees' -> forest_minimize

        args: 
            (1) optimization_method (type:str) - one of the optimization method to employ

        return: 
            (1) Validation loss
        """

        @use_named_args(dimensions = self.dimensions)
        def objective(learning_rate, num_dense_layers, num_dense_nodes, batch_size):
            """
            Objective function to minimize
            """

            # Create ANN model
            model = self.create_model(
                                    learning_rate = learning_rate,
                                    num_dense_layers = num_dense_layers,
                                    num_dense_nodes = num_dense_nodes,
                                    batch_size = batch_size 
                                    )


            # Fit to data
            if self.problem == 'classification':
                y_train = tf.keras.utils.to_categorical(
                                    self.sample_train[self.target], 
                                    num_classes = self.target_classes
                                    )

                blackbox = model.fit(
                                    x = self.sample_train[self.RFE_features],
                                    y = y_train,
                                    epochs = 100, 
                                    batch_size = batch_size*8,
                                    validation_split = 0.2,
                                    shuffle = True,
                                    callbacks = self.callback
                                    )


            elif self.problem == 'regression':
                blackbox = model.fit(
                                    x = self.sample_train[self.RFE_features],
                                    y = self.sample_train[self.target],
                                    epochs = 100,
                                    batch_size = batch_size*8,
                                    validation_split = 0.2,
                                    shuffle = True,
                                    callbacks = self.callback
                                    )
        

            # Return the validation accuracy for the last epoch.
            val_loss = blackbox.history['val_loss'][-1]
            
            print()
            print("loss: {0:.2f}".format(val_loss))
            print()


            # Delete the Keras model with these hyper-parameters from memory.
            del model

            
            # Clear the Keras session
            K.clear_session()
            tf.compat.v1.reset_default_graph()

            
            return val_loss


        # Apply optimization technique
        self.optimization_method = optimization_method

        if self.optimization_method == 'random_search':
            opt_method = dummy_minimize

        elif self.optimization_method == 'bayesian':
            opt_method = gp_minimize

        elif self.optimization_method == 'gradient_bossted_trees':
            opt_method = gbrt_minimize

        elif self.optimization_method == 'decision_trees':
            opt_method = forest_minimize


        self.opt = opt_method(
                            func = objective,
                            dimensions = self.dimensions,
                            n_calls = 100, 
                            n_initial_points = 10,
                            random_state = 42,
                            #n_jobs = -1,
                            noise = 'gaussian',
                            x0 = self.default_parameters,
                            verbose = 1
                            )

        self.values = list()

        print('\n', '*** Optimal hyperparameters *** ')

        for i in range(0, len(self.opt.x)): 
            print('{}: {}'.format(self.hyperparameters[i], self.opt.x[i]))
            self.values.append(self.opt.x[i])

        dump(opt_method, os.path.join(self.path_to_save, r'optimization_ann_data.pkl'))

        print('Optimisation data saved as: "optimization_ann_data.pkl"')



    def convergence_plot(self):
        """
        Plot convergence plot of the optimization
        """

        # Setting up the figure
        fig, ax = plt.subplots(figsize = (8,8))

        fontsize = 16

        plot = plot_convergence((str(self.optimization_method), self.opt))

        plot.legend(loc = "best", prop = {'size': fontsize}, numpoints = 1)
        ax.grid(b = None)
        ax.set_title(' ', fontsize = 18)
        ax.set_xlabel('Number of iterations', fontsize = fontsize)
        ax.set_ylabel('Objective minimum', fontsize = fontsize) 
        ax.tick_params(axis = 'both', which = 'major', labelsize = fontsize, direction = 'in')


        #final_figure
        fig.savefig(os.path.join(self.path_to_save, r'Optimisation_result.png'), dpi = 300, bbox_inches = "tight")

        print('Figure saved as: "Optimisation_result.png"')



    def objective_plot(self):
        """
        Plot objective and corresponding evaluation plots
        """
        _ = plot_objective(self.opt, n_points = 10) 
        _ = plot_evaluations(self.opt)



    def train_model(self):
        """
        Train model with optimal hyperparameters identified during optimiation
        """

        # Clear sessionn
        K.clear_session()
        tf.compat.v1.reset_default_graph()   


        # Naming convention to use for saved weights
        checkpoint_name = 'Weights-{epoch:03d}--{loss:.5f}.hdf5' 


        # Define callback 
        if self.problem == 'regression':
            callback2 = [EarlyStopping(
                                    monitor = 'loss', 
                                    patience = 5, 
                                    min_delta = 0, 
                                    verbose = 1, 
                                    mode = 'min'
                                    ),

                        ModelCheckpoint(
                                    checkpoint_name, 
                                    monitor = 'loss', 
                                    verbose = 1, 
                                    save_best_only = True, 
                                    mode = 'min'
                                    )] 


        elif self.problem == 'classification':
            callback2 = [EarlyStopping(
                                    monitor = 'loss', 
                                    patience = 5, 
                                    min_delta = 0, 
                                    verbose = 1, 
                                    mode = 'max'
                                    ),

                        ModelCheckpoint(
                                    checkpoint_name, 
                                    monitor = 'loss', 
                                    verbose = 1, 
                                    save_best_only = True, 
                                    mode = 'min'
                                    )] 


        # Optimal parameters
        learning_rate = self.opt.x[0]
        num_dense_layers = self.opt.x[1]
        num_dense_nodes = self.opt.x[2]
        batch_size = self.opt.x[3]

        # Create ANN model using optimal hyperparameters
        self.model = self.create_model(
                                    learning_rate = learning_rate,
                                    num_dense_layers = num_dense_layers,
                                    num_dense_nodes = num_dense_nodes,
                                    batch_size = batch_size
                                    )

        # Fit to training data
        if self.problem == 'classification':
            self.y_train = tf.keras.utils.to_categorical(
                                self.sample_train[self.target], 
                                num_classes = self.target_classes
                                )

            self.y_test= tf.keras.utils.to_categorical(
                                self.sample_test[self.target], 
                                num_classes = self.target_classes
                                )

            self.model.fit(
                        x = self.sample_train[self.RFE_features], 
                        y = self.y_train, 
                        epochs = 10, #300, 
                        batch_size = batch_size*8, 
                        callbacks = callback2,
                        verbose = 1
                        )

            # Evaluate using test set
            self.model.evaluate(
                                self.sample_test[self.RFE_features], 
                                self.y_test
                                )


        # Fit to training data
        elif self.problem == 'regression':
            self.model.fit(
                        x = self.sample_train[self.RFE_features], 
                        y = self.sample_train[self.target], 
                        epochs = 10, #300, 
                        batch_size = batch_size*8, 
                        callbacks = callback2,
                        verbose = 1
                        )

            # Evaluate using test set
            self.model.evaluate(
                                self.sample_test[self.RFE_features], 
                                self.sample_test[self.target]
                                )


        return self.model.summary()



    def choose_weight(self, path_to_weight):
        """
        Choose the weight for the ANN model

        args: 
            (1) path_to_weight (type:str) - path to the latest weight file e.g. 'Weights-083--13.58393.hdf5'

        return:
            (1) model with selected weight (note: model != self.model)
        """

        # Optimal parameters
        learning_rate = self.opt.x[0]
        num_dense_layers = self.opt.x[1]
        num_dense_nodes = self.opt.x[2]
        batch_size = self.opt.x[3]


        # Load wights file of the best model
        wights_file = path_to_weight 
        self.model.load_weights(wights_file) 


        # Define learning rate schedule
        lr_schedule = ExponentialDecay(
                                        initial_learning_rate = learning_rate,
                                        decay_steps = 100000,
                                        decay_rate = 0.96
                                        )

        # Define optimizer
        optimizer = Adam(learning_rate = lr_schedule)


        # Compile model
        if self.problem == 'regression':
            self.model.compile(
                            optimizer = optimizer,
                            loss = 'mean_squared_error',
                            metrics = [
                                        metrics.MeanSquaredError(), 
                                        metrics.RootMeanSquaredError(), 
                                        metrics.MeanAbsoluteError()
                                        ])


        elif self.problem == 'classification':
            if self.multiclass == True:
                self.model.compile(
                            loss = 'categorical_crossentropy', 
                            optimizer = optimizer, 
                            metrics = [
                                        metrics.Accuracy(name = 'accuracy'),
                                        metrics.Precision(name = 'precision'),
                                        metrics.Recall(name = 'recall')
                                        ])


            elif self.multiclass == False:
                self.model.compile(
                        loss = 'binary_crossentropy', 
                        optimizer = optimizer, 
                        metrics = [
                                    metrics.Accuracy(name = 'accuracy'),
                                    metrics.Precision(name = 'precision'),
                                    metrics.Recall(name = 'recall')
                                    ])


        return self.model



    def regression_plot(self, X, Y, min_value, max_value):
        """
        Show regression results
        This function is recalled using the 'evaluate()' function

        args: 
            (1) X (type:list) - true/observed target values
            (2) Y (type:list) - predicted target values
            (3) min_value (type:int) - min value to plot i.e. lower limit
            (4) max_value (type:int) - max value to plot i.e. upper limit

        return: 
            (1) stats and figure of regression analysis
        """

        # Figure
        plt.figure(figsize = (8, 8)) 


        # Predicted vs Actual
        plt.plot(X, Y, 'o', markersize = 5, color = 'black', alpha = 0.15)


        # line of best fit
        no_ticks = max_value
        linear_fit = np.linspace(0, no_ticks - 5, no_ticks)
        plt.plot(linear_fit, linear_fit*self.stats_results.params[1] + self.stats_results.params[0], '-', color = 'tab:blue') 


        # Ideal y=x 
        y = x = np.linspace(0, no_ticks - 5, no_ticks)
        plt.plot(x, y, '--', color='red', alpha=0.8) 


        fontsize = 18
        plt.xlim([min_value, no_ticks])
        plt.ylim([min_value, no_ticks])
        plt.xlabel('True target value', fontsize=fontsize)
        plt.ylabel('Predicted target value', fontsize=fontsize)
        plt.tick_params(axis='both', which='both', labelsize=fontsize, direction="in")
        plt.rcdefaults()

        print('Linear fit has: ')
        print('m = ', self.stats_results.params[1])
        print('c = ', self.stats_results.params[0], '\n')

        plt.savefig(os.path.join(self.path_to_save, r'regression_plot.png'), dpi = 300, bbox_inches = "tight")
        plt.show()

        print('Figure saved as: "regression_plot.png"')



    def confusion_matrix(self, target_names):
        """
        Generate confusion matrix plot 
        This function is recalled using the 'evaluate()' function

        args: 
            (1) target_names (type:list) = list of target classes e.g. ['AFM', 'FM', 'NM']
        
        return: 
            (1) conusion matrix plot
        """

        # Pretty confusion matrix 
        disp = plot_confusion_matrix(self.model, 
                                     self.sample_test[self.RFE_features], 
                                     self.sample_test[self.target],
                                     display_labels = np.array(target_names, dtype = '<U10'),
                                     cmap = plt.cm.Blues,
                                     normalize = None)
        

        fontsize = 13
        plt.tick_params(axis = 'both', which = 'major', labelsize = fontsize, direction = 'in')
        plt.savefig(os.path.join(self.path_to_save, r'Confusion_matrix.png'), dpi = 300, bbox_inches = 'tight')
        plt.show()

        print('Figure saved as: "Confusion_matrix.png"')



    def evaluate(self, strategy, *args, **kwargs):
        """
        Evaluate the ML model using out-of-sample test set

        args: 
            (1) strategy (type:str); averaging method e.g. 'micro', 'macro', 'weighted'
            (2*) target_names (type:list) = list of target classes e.g. ['AFM', 'FM', 'NM']

        return: 
            (1) stats and plots of result
        """

        if self.problem == 'classification':

            target_names = kwargs.get('target_names')


            # Apply model onto test data
            self.y_test = self.sample_test[self.target]


            # In tf, model.predict() returns probability
            self.y_pred =  self.model.predict(self.sample_test[self.RFE_features])
            self.y_pred_2 = np.argmax(self.y_pred, axis = 1)


            # Evaluate metric scores
            print('1. The F-1 score of the model {}\n'.format(f1_score(self.y_test.ravel(), self.y_pred_2, average = strategy)))
            print('2. The recall score of the model {}\n'.format(recall_score(self.y_test.ravel(), self.y_pred_2, average = strategy)))
            print('3. Classification report \n {} \n'.format(classification_report(self.y_test.ravel(), self.y_pred_2, target_names = target_names)))
            print('4. Classification report \n {} \n'.format(multilabel_confusion_matrix(self.y_test.ravel(), self.y_pred_2)))
            print('5. Confusion matrix \n {} \n'.format(confusion_matrix(self.y_test.ravel(), self.y_pred_2)))
            print('6. Accuracy score \n {} \n'.format(accuracy_score(self.y_test.ravel(), self.y_pred_2)))
            print('7. Balanced accuracy score \n {} \n'.format(balanced_accuracy_score(self.y_test.ravel(), self.y_pred_2)))


            # Evaluate matthews correlation coef
            y_test_2 = label_binarize(self.y_test, classes=[i for i in range(self.target_classes)])


            # Convert each row to 1 and 0 based on prob
            all_scores = self.y_pred
            all_scores_2 = np.zeros_like(all_scores)
            all_scores_2[np.arange(len(all_scores)), all_scores.argmax(1)] = 1


            print('8. Matthews corrcoef of Class: ')

            m_corr = list()

            for i in range(self.target_classes):
                corr = matthews_corrcoef(y_test_2[:, i], all_scores_2[:, i])
                m_corr.append(corr)
                print(str(target_names[i]) + ': ', corr)


            print('9. Matthews macro corrcoef \n {} \n'.format(sum(m_corr)/3)) 


        elif self.problem == 'regression':

            adjusted = kwargs.get('adjusted')
            min_value = kwargs.get('min_value')
            max_value = kwargs.get('max_value')


            # Apply model onto test data
            self.y_test = self.sample_test[self.target]
            self.y_pred = self.model.predict(self.sample_test[self.RFE_features]).ravel()
            self.id_index = self.sample_test.index.tolist()


            df_pred = pd.DataFrame(
                                    {'task_id': self.id_index, 
                                    str(self.target): self.y_test, 
                                    'pred_target': self.y_pred
                                    })


            # Create a column to eliminate negative values
            df_pred['adjusted_pred_target'] = df_pred['pred_target']
            df_pred['adjusted_pred_target'] = df_pred['adjusted_pred_target'].apply(lambda x: 0 if x < 0 else x)

            # Define X and Y variable
            X = df_pred[self.target]

            if adjusted == True:
                # use adjusted values if negative value of the feature does not make sense
                Y = df_pred['adjusted_pred_target']

            else:
                Y = df_pred['pred_target'] 


            # Stats
            self.stats_results = sm.OLS(Y,sm.add_constant(X)).fit()

            print(self.stats_results.summary())

            print('MAE: ', mean_absolute_error(X, Y))
            print('MSE: ', mean_squared_error(X, Y))
            print('RMSE: ', mean_squared_error(X, Y, squared = False))
            print('R-squared: ', r2_score(X, Y))
            print('Max error: ', max_error(X, Y))
            print('Explained_variance_score: ', explained_variance_score(X, Y, multioutput = 'variance_weighted'))


            # Plot figure
            self.regression_plot(X, Y, min_value, max_value)



    def ROC(self, overall_performance, *args, **kwargs):
        """
        Generate ROC plot for the classification problem

        args: 
            (1) overall_performance (type:bool); whether to plot the overall average, where strategy determines the method of averaging
            (2*) strategy (type:str); averaging method e.g. 'micro', 'macro', 'weighted'

        return: 
            (1) figure of ROC 
        """

        strategy = kwargs.get('strategy')

        self.y_test = self.sample_test[self.target]
        self.y_pred =  self.model.predict(self.sample_test[self.RFE_features])
        self.y_pred_2 = np.argmax(self.y_pred, axis = 1)


        # Compute ROC curve and ROC area for each class
        self.fpr = dict()
        self.tpr = dict()

        n_classes = self.y_pred.shape[1]
        roc_auc = dict()

        self.y_test_2 = label_binarize(self.y_test, classes = list(range(n_classes)))


        #################### Micro
        for i in range(n_classes):
            self.fpr[i], self.tpr[i], _ = roc_curve(self.y_test_2[:, i], self.y_pred[:, i])
            roc_auc[i] = auc(self.fpr[i], self.tpr[i])

        # Compute micro-average ROC curve and ROC area
        self.fpr["micro"], self.tpr["micro"], _ = roc_curve(self.y_test_2.ravel(), self.y_pred.ravel())
        roc_auc["micro"] = auc(self.fpr["micro"], self.tpr["micro"])


        #################### Macro
        # First aggregate all false positive rates
        all_fpr = np.unique(np.concatenate([self.fpr[i] for i in range(n_classes)]))

        # Then interpolate all ROC curves at this points
        mean_tpr = np.zeros_like(all_fpr)
        for i in range(n_classes):
            mean_tpr += np.interp(all_fpr, self.fpr[i], self.tpr[i])


        # Finally average it and compute AUC
        mean_tpr /= n_classes

        self.fpr["macro"] = all_fpr
        self.tpr["macro"] = mean_tpr
        roc_auc["macro"] = auc(self.fpr["macro"], self.tpr["macro"])


        # Plot all ROC curves
        plt.figure(figsize=(8,8))

        if overall_performance == True:
            if strategy == 'micro':
                plt.plot(
                        self.fpr["micro"], self.tpr["micro"],
                        label = 'micro-average ROC (AUC = {0:0.3f})'
                               ''.format(roc_auc["micro"]),
                        color = 'tab:green', 
                        linestyle = '-', 
                        linewidth = 4
                        )


            if strategy == 'macro':
                plt.plot(
                        self.fpr["macro"], self.tpr["macro"],
                        label = 'macro-average ROC (AUC = {0:0.3f})'
                            ''.format(roc_auc["macro"]),
                        color = 'tab:blue', 
                        linestyle = '-', 
                        inewidth = 4
                        )


        if overall_performance == False:
            # Individual class
            lw = 2
            colors = cycle(['aqua', 'darkorange', 'cornflowerblue'])

            for i, color in zip(range(n_classes), colors):
                plt.plot(
                        self.fpr[i], 
                        self.tpr[i], 
                        color = color, 
                        lw = lw,
                        label = 'ROC curve of class {0} (AUC = {1:0.3f})'
                            ''.format(i, roc_auc[i])
                        )


        # Plot curves
        fontsize = 18
        lw = 2

        plt.plot([0, 1], [0, 1], 'k--', lw = lw)
        plt.xlim([0.0, 1])
        plt.ylim([0.0, 1.01])
        plt.xlabel('False Positive Rate', fontsize = fontsize)
        plt.ylabel('True Positive Rate', fontsize = fontsize)
        plt.tick_params(axis = 'both', which = 'major', labelsize = fontsize)
        plt.legend(loc = "lower right", fontsize = fontsize, framealpha = 1)


        #final_figure
        plt.savefig(os.path.join(self.path_to_save, r'Receiver_operating_characteristic_curve.png'), dpi = 300, bbox_inches = "tight")
        plt.show()

        print('Figure saved as: "Receiver_operating_characteristic_curve.png"')



    def DET(self, strategy):
        """
        Generate DET plot for the classification problem

        args: 
            (1) strategy (type:str); averaging method e.g. 'micro', 'macro', 'weighted'

        return: 
            (1) figure of DET curve
        """

        # Detection Error Trade-off Curve
        fnr_macro = 1 - self.tpr['macro']
        fnr_micro = 1 - self.tpr['micro']


        # Plot curves
        fontsize = 18
        linewidth = 2


        plt.figure(figsize = (8,8))

        if strategy == 'macro':
            plt.plot(
                    fnr_macro, 
                    self.fpr['macro'] ,
                    label = 'macro-average ERR',
                    color = 'tab:blue', 
                    linestyle = '-',
                    linewidth = linewidth
                    )


        if strategy == 'micro':
            plt.plot(
                    fnr_micro, 
                    self.fpr['micro'] ,
                    label = 'micro-average ERR ',
                    color = 'tab:green', 
                    linestyle = '-',
                    linewidth = linewidth
                    )


        lw = 2
        plt.plot([0, 1], [0, 1], 'k--', lw = lw)
        plt.xlim([0.0, 1])
        plt.ylim([0.0, 1])
        plt.xlabel('False Negative Rate', fontsize = fontsize)
        plt.ylabel('False Positive Rate', fontsize = fontsize)
        plt.tick_params(axis = 'both', which = 'major', labelsize = fontsize, direction = 'in')
        plt.legend(loc = "upper right", fontsize=fontsize, framealpha=1)


        #final_figure
        plt.savefig(os.path.join(self.path_to_save, r'detection_error_tradeoff_curves_v1.png'), dpi = 300, bbox_inches = "tight")
        plt.show()

        print('Figure saved as: "detection_error_tradeoff_curves_v1.png"')

    

    def PR(self):
        """
        Generate PR curve for the classification problem
        """

        self.y_test = self.sample_test[self.target]
        self.y_pred =  self.model.predict(self.sample_test[self.RFE_features])
        self.y_pred_2 = np.argmax(self.y_pred, axis = 1)


        # For each class
        n_classes = self.y_pred.shape[1]
        precision = dict()
        recall = dict()
        average_precision = dict()
        thresholds = dict()
        
        self.y_test_2 = label_binarize(self.y_test, classes = list(range(n_classes))) 
            

        # For each class / for the top classifier
        for i in range(n_classes):
            precision[i], recall[i], _ = precision_recall_curve(self.y_test_2[:, i], self.y_pred[:, i])
            average_precision[i] = average_precision_score(self.y_test_2[:, i], self.y_pred[:, i])
                

        precision["micro"], recall["micro"], thresholds['micro'] = precision_recall_curve(self.y_test_2.ravel(),self.y_pred.ravel())
            

        average_precision["micro"] = average_precision_score(self.y_test_2, self.y_pred, average = "micro")
        average_precision["weighted"] = average_precision_score(self.y_test_2, self.y_pred, average = "weighted")
        average_precision["macro"] = average_precision_score(self.y_test_2, self.y_pred, average = "macro")
    

        print('Average precision score, micro-averaged over all classes: {0:0.3f}'
                      .format(average_precision["micro"]))

        print('Average precision score, macro-averaged over all classes: {0:0.3f}'
                      .format(average_precision["macro"]))

        print('Average precision score, weighted-averaged over all classes: {0:0.3f}'
                      .format(average_precision["weighted"]))

        #print('PR_AUC_micro: ', auc(recall["micro"], precision["micro"]))


        # Plot figure
        plt.figure(figsize = (8,8))

        fontsize = 18

        plt.step(
                recall['micro'], precision['micro'], 
                where = 'post', 
                lw = 2, 
                color = 'tab:blue', 
                label = 'Micro-averaged PR (AP = 0.995)'
                )


        labelsize = 18

        plt.xlabel('Recall',fontsize = fontsize)
        plt.ylabel('Precision',fontsize = fontsize)
        plt.ylim([0.0, 1.05])
        plt.xlim([0.0, 1.0])
        plt.tick_params(axis = 'both', which = 'major', labelsize = labelsize, direction = 'in')
        plt.legend(fontsize = fontsize, loc = "lower left", framealpha = 1.0)


        f_scores = np.linspace(0.2, 0.8, num = 4)
        lines = []
        labels = []

        n = 0
        for f_score in f_scores:
            x = np.linspace(0.001, 1.0)
            y = f_score * x / (2 * x - f_score)
            l, = plt.plot(x[y >= 0], y[y >= 0], color = 'gray', alpha = 0.2)


        # Location of the annotation
        x0 = [0.13, 0.26, 0.43, 0.67]
        y0 = [0.2, 0.4, 0.6, 0.8]
        n = 0

        fontsize2 = 14

        while n < len(x0):
            if n < 0:
                plt.annotate('F1 =' + str(y0[n]), xy = (x0[n], 0.99 + 0.02),fontsize = fontsize2)
            else:
                plt.annotate('F1 =' + str(y0[n]), xy = (x0[n], 0.99 + 0.02),fontsize = fontsize2)
            n = n + 1 


        #Save figure
        plt.savefig(os.path.join(self.path_to_save, r'precision_recall.png'), dpi = 300, bbox_inches = "tight")
        plt.show()

        print('Figure saved as: "precision_recall.png"')

