import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.neural_network import MLPClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import f1_score
from imblearn.metrics import geometric_mean_score
from sklearn.feature_selection import r_regression, f_regression
from sklearn.datasets import make_classification
from collections import Counter
from scipy.stats import spearmanr, ConstantInputWarning
from imbalance_measures import measure_imbalance, status_of_imbalance, multiclass_imbalance_ratio
import seaborn as sns
import sys
import warnings
from sklearn.exceptions import ConvergenceWarning
from pycol_complexity.complexity import Complexity
from pymfe.mfe import MFE
import math
pd.options.mode.chained_assignment = None

def make_DS_descriptive_table(all_ds, col_label='label', only_class_freq=False):
    """
    Makes a descriptive table of the datasets contained in the data/ folder and saves it as a .csv
    in the table/ folder.
    @all_ds: list of Strings, list of the names of all datasets that should be in the table 
    @col_label: (optional) String, default='label', name of the labels column in the datasets 
                [WARNING]: the name of the labels column should be the same in all the datasets!
    @only_class_freq: (optional) boolean, default=False, 
                      if False, the table contains: dataset name, number of classes, number of features, total number
                                                    of points, mean and sd of class frequencies, results of all 8
                                                    imbalance metrics (IR, adjIR, ID, LRID, IF, C1, C2, SIMBA)
                      if True, the table contains:  dataset name, number of classes, number of minority classes,
                                                    list of class frequencies
    """
    # Initialization of Dataframe with correct columns
    if not only_class_freq:
        descriptive_table = pd.DataFrame(columns=["Dataset", "#c", "#ft", "Size", "Class frequencies (%) mean ± sd", 
                                                  "IR", "Adj-IR", "C1", "C2", "ID", "LRID", "IF", "SIMBA"])
    else:
        descriptive_table = pd.DataFrame(columns=["Dataset", "#c", "#min", "Class frequencies (%)"])

    for ds in all_ds: # for each dataset
        path = "./data/" + ds + "/"
        df = pd.read_csv(path + "df_allfeatures.csv") 
        labels = pd.read_csv(path + "labels.csv")
        name = ds
        print("NAME: " + name)
        n_class = len(labels[col_label].unique())
        n_features = len(df.columns) - 1 # - 1 to account for the column of id
        n_samples = len(labels)
        n_class_freq = Counter(labels[col_label])
        class_freq_percentage = []
        n_min_class = 0
        for i in range(len(n_class_freq.keys())):
            percentage = round((n_class_freq[i]/n_samples)*100,2) # class freq is rounded to 2 decimals
            if percentage < 100/n_class:
                n_min_class += 1 # a class is considered a minority class if its freq is < 100/n_class
            class_freq_percentage.append(percentage)
        df = df.drop(columns=['id'])
        df[col_label] = labels[col_label]
        if not only_class_freq:
            # imbalance metrics calculation
            IR, adjIR, C1, C2, ID, LRID, IF, SIMBA = measure_imbalance(df, col_label, rounded=True)
            mean_cf = round(np.mean(class_freq_percentage),2)
            sd_cf = round(np.std(class_freq_percentage),2)
        if not only_class_freq:
            row = [name, n_class, n_features, n_samples, str(mean_cf) + "±" + str(sd_cf), IR, adjIR,  C1, C2, 
                   ID, LRID, IF, SIMBA]
        else:
            row = [name, n_class, n_min_class, class_freq_percentage]
        descriptive_table.loc[len(descriptive_table.index)] = row # add the row corresponding to the current dataset
    print(descriptive_table)
    # saving the table in table/ folder
    if not only_class_freq:
        descriptive_table.to_csv("./table/ds_descriptivetable.csv", index=False)
    else:
        descriptive_table.to_csv("./table/ds_descriptivetable_classfreq.csv", index=False)

def get_eval_score(X_train, y_train, X_test, y_test, classif, i, eval_method):
    """
    Performs classification and returns the macro average f1-score (float) or the geometric mean (float)
    @X_train, X_test: DataFrames, training and testing datasets (does not contain labels)
    @y_train, y_test: Series, it corresponds to the labels column from the dataset,  
                      training and testing labels (i.e. ground truth)
    @classif: String, chosen classifier, can take the following values: "SVM", "LDA", "kNN", "RF", and "MLP"
              [WARNING!] any other value will lead to an error and the program will be stopped.
    @i: integer, random state value. 
    @eval_method: String, defines the classification evaluation metric: either "f1" for f1-score or gmean
    """
    # train
    if classif == "SVM":
        clf = make_pipeline(StandardScaler(), SVC(probability=True))
    elif classif == "LDA":
        clf = make_pipeline(StandardScaler(), LinearDiscriminantAnalysis())
    elif classif =="kNN":
        neighb = int(np.sqrt(len(X_train)))
        clf = make_pipeline(StandardScaler(), KNeighborsClassifier(n_neighbors=neighb))
    elif classif == "RF":
        clf = make_pipeline(StandardScaler(), RandomForestClassifier(random_state=i))
    elif classif == "MLP":
        warnings.simplefilter('ignore', category=ConvergenceWarning)
        clf = make_pipeline(StandardScaler(), MLPClassifier(random_state=i))
    else:
        print("[ERROR] Unknown specified classifier.")
        sys.exit()
    clf.fit(X_train, y_train)
    # test
    y_pred = clf.predict(X_test)
    if eval_method == "f1":
        eval_score = f1_score(y_test, y_pred, average='macro')
    else:
        eval_score = geometric_mean_score(y_test, y_pred, average='macro')
    return eval_score

def get_metrics_eval_score(data, labels, classif="SVM", onlyMetrics=False, metric="all", eval_method="f1"):
    """
    Returns the avg of all imbalance metrics (IR, adjIR, C1, C2, ID, LRID, IF, and SIMBA) or of SIMBA only and the avg f1-score/gmean 
    score obtained with 5-CV.
    @data: pd.DataFrame, containing all features (not labels)
    @labels: pd.DataFrame, containing labels
    @classif: (optional) String, default="SVM", chosen classifier, can take the following values: "SVM", "LDA", 
              "kNN", "RF", and "MLP". Will only be used if onlyMetrics=False. 
              [WARNING!] any other value will lead to an error and the program will be stopped.
    @onlyMetrics: (optional) boolean, default=False, if True avg f1-score/gmean will not be calculated and put to 0. 
    @metric: String, if "all", all metrics returned, otherwise, only SIMBA is returned.
    @eval_method: String, defines the classification evaluation metric: either "f1" for f1-score or gmean
    """
    l_evalscores = []
    l_IR = []
    l_adjIR = []
    l_C1 = []
    l_C2 = []
    l_ID = []
    l_LRID = []
    l_IF = []
    l_SIMBA = []
    X = data
    y = labels
    kf = StratifiedKFold(n_splits=5) # 5-fold CV
    for i, (train_index, test_index) in enumerate(kf.split(X, y)):
        # print(f"Fold {i + 1}.")
        X_test = X.iloc[test_index]
        y_test = y.iloc[test_index]
        X_train = X.iloc[train_index]
        y_train = y.iloc[train_index]
        # metric calculations
        col_label = str(y_train.name)
        df = X_train.copy()
        df[col_label] = y_train
        if metric == "all":
            IR, adjIR, C1, C2, ID, LRID, IF, SIMBA = measure_imbalance(df, col_label)
            l_IR.append(IR)
            l_adjIR.append(adjIR)
            l_C1.append(C1)
            l_C2.append(C2)
            l_ID.append(ID)
            l_LRID.append(LRID)
            l_IF.append(IF)
            l_SIMBA.append(SIMBA)
        else:
            l_SIMBA.append(status_of_imbalance(df, col_label))
        if not onlyMetrics:
            eval_score = get_eval_score(X_train, y_train, X_test, y_test, classif, i, eval_method=eval_method)
        else:
            eval_score = 0
        l_evalscores.append(eval_score)
    # Averaging the results of the 5 runs
    if metric == "all":
        avg_IR = np.mean(l_IR)
        avg_adjIR = np.mean(l_adjIR)
        avg_C1 = np.mean(l_C1)
        avg_C2 = np.mean(l_C2)
        avg_ID = np.mean(l_ID)
        avg_LRID = np.mean(l_LRID)
        avg_IF = np.mean(l_IF)
    avg_SIMBA = np.mean(l_SIMBA)
    avg_eval_score = np.mean(l_evalscores)
    if metric == "all":
        return avg_IR, avg_adjIR, avg_C1, avg_C2, avg_ID, avg_LRID, avg_IF, avg_SIMBA, avg_eval_score
    else:
        return avg_SIMBA, avg_eval_score

def get_SIMBA_ablation_eval_score(data, labels, classif="SVM", onlyMetrics=False, eval_method="f1"):
    """
    Returns the avg of all versions of SIMBA (whole, nonorm, noftimp, noftred) with 5-CV for ablation study.
    @data: pd.DataFrame, containing all features (not labels)
    @labels: pd.DataFrame, containing labels
    @classif: (optional) String, default="SVM", chosen classifier, can take the following values: "SVM", "LDA", 
              "kNN", "RF", and "MLP". Will only be used if onlyMetrics=False. 
              [WARNING!] any other value will lead to an error and the program will be stopped.
    @onlyMetrics: (optional) boolean, default=False, if True avg f1-score/gmean will not be calculated and put to 0. 
    @eval_method: String, defines the classification evaluation metric: either "f1" for f1-score or gmean
    """
    l_evalscores = []
    l_SIMBA_whole = []
    l_SIMBA_nonorm = []
    l_SIMBA_noftimp = []
    l_SIMBA_noftred = []
    X = data
    y = labels
    kf = StratifiedKFold(n_splits=5) # 5-fold CV
    for i, (train_index, test_index) in enumerate(kf.split(X, y)):
        # print(f"Fold {i + 1}.")
        X_test = X.iloc[test_index]
        y_test = y.iloc[test_index]
        X_train = X.iloc[train_index]
        y_train = y.iloc[train_index]
        # metric calculations
        col_label = str(y_train.name)
        df = X_train.copy()
        df[col_label] = y_train
        l_SIMBA_whole.append(status_of_imbalance(df, col_label))
        l_SIMBA_nonorm.append(status_of_imbalance(df, col_label, version='nonorm'))
        l_SIMBA_noftimp.append(status_of_imbalance(df, col_label, version='noftimp'))
        l_SIMBA_noftred.append(status_of_imbalance(df, col_label, version='noftred'))
        if not onlyMetrics:
            eval_score = get_eval_score(X_train, y_train, X_test, y_test, classif, i, eval_method=eval_method)
        else:
            eval_score = 0
        l_evalscores.append(eval_score)
    # Averaging the results of the 5 runs
    avg_SIMBA_whole = np.mean(l_SIMBA_whole)
    avg_SIMBA_nonorm = np.mean(l_SIMBA_nonorm)
    avg_SIMBA_noftimp = np.mean(l_SIMBA_noftimp)
    avg_SIMBA_noftred = np.mean(l_SIMBA_noftred)
    avg_eval_score = np.mean(l_evalscores)
    return avg_SIMBA_whole, avg_SIMBA_nonorm, avg_SIMBA_noftimp, avg_SIMBA_noftred, avg_eval_score

def add_correlation(evalscores, all_metric_values, avg_corr_metric, l_corr_metric):
    # Function used to solve an issue with the r_regression function. 
    # Sometimes r_regression returns a number outside the range [-1,+1] when PCC==0; thus this is solved as follows
    if abs(r_regression(evalscores, all_metric_values)) <= 1:
        avg_corr_metric += r_regression(evalscores, all_metric_values)
        l_corr_metric.append(r_regression(evalscores, all_metric_values))
    else:
        l_corr_metric.append(0)
    return avg_corr_metric, l_corr_metric

def run_synthdata(N = 1, data_distrib = [100, 5, 0], feature_distrib = [2, 0, 8], tot_ft = 10, i_range=[5, 101, 1], 
              title='Synthetic datasets', corr=True, varying="data", eval_method="f1"):
    """
    Creates and runs the tests on synthetic datasets. Either computes metrics and f1-scores/gmeans for each dataset, and 
    plots trend of the metric, PCC and SRCC between metrics scores and f1-scores/gmeans (corr=True). Or, computes only 
    the metrics and return them, but no plot is done (corr=False). The dataset generated can vary on data
    distribution or on the feature distribution (informative, redundant, not informative).
    @N: integer, number of random runs. For each run, a different random state is picked. 
    @data_distrib: list of integers size 3, we use datasets with 3 classes. data_distrib gives the number of data
                   points in each class. If one element of the list == 0, this element is the varying one (it will
                   take all values in i_range). 
    @feature_distrib: list of integers size 3, this represents [#informative ft, #redundant ft, #not informative ft],
                      if varying is not on data but on ft, one of these might be equal to "-1" and will therefore be
                      the varying element. tot_ft can also vary according to this.
    @tot_ft: integer, total number of features. It is equal to sum(feature_distib) when data distribution is varying.
             If ft is varying, then it is possible that it is not fixed (i.e. == -1), then it will vary with i.   
    @i_range: list of integers, size 3, i is the varying factor, whether for data or feature distribution. i will 
              take any values in the specified range [min (incl.), max (excl.), step].
    @title: String, title of the plots, only used if corr==True.
    @corr: boolean, if True, metrics and f1-scores/gmeans are computed for each dataset for each i, and 3 plots are made:
              trend of the metrics when i is varying, PCC and SRCC between metrics scores and f1-scores/gmeans. If False, the
              f1-scores/gmeans are not computed, only the metrics scores, which are returned. 
    @varying: String, defines what is varying between "data" or features. If == "data", data distribution is varying.
              If != "data", it is considered that the features are varying.
    @eval_method: String, defines the classification evaluation metric: either "f1" for f1-score or gmean
    """
    print('In process: ' + title)

    avg_all_IR = []
    avg_all_adjIR = []
    avg_all_C1 = []
    avg_all_C2 = []
    avg_all_ID = []
    avg_all_LRID = []
    avg_all_IF = []
    avg_all_SIMBA = []
    avg_eval_scores = []
    avg_index = []
    if corr:
        avg_corr_IR = 0
        avg_corr_adjIR = 0
        avg_corr_C1 = 0
        avg_corr_C2 = 0
        avg_corr_ID = 0
        avg_corr_LRID = 0
        avg_corr_IF = 0
        avg_corr_SIMBA = 0
        l_corr_IR = []
        l_corr_adjIR = []
        l_corr_C1 = []
        l_corr_C2 = []
        l_corr_ID = []
        l_corr_LRID = []
        l_corr_IF = []
        l_corr_SIMBA = []
        avg_corr_IR_srcc = 0
        avg_corr_adjIR_srcc = 0
        avg_corr_C1_srcc = 0
        avg_corr_C2_srcc = 0
        avg_corr_ID_srcc = 0
        avg_corr_LRID_srcc = 0
        avg_corr_IF_srcc = 0
        avg_corr_SIMBA_srcc = 0
        l_corr_IR_srcc = []
        l_corr_adjIR_srcc = []
        l_corr_C1_srcc = []
        l_corr_C2_srcc = []
        l_corr_ID_srcc = []
        l_corr_LRID_srcc = []
        l_corr_IF_srcc = []
        l_corr_SIMBA_srcc = []

    for k in range(N): # each random run/trial
        print("Trial: " + str(k+1) + "/" + str(N))
        all_IR = []
        all_adjIR = []
        all_C1 = []
        all_C2 = []
        all_ID = []
        all_LRID = []
        all_IF = []
        all_SIMBA = []
        eval_scores = []
        index = []
        for i in range(i_range[0], i_range[1], i_range[2]): # iterating over i values
            if varying == "data":
                # Varying data distribution
                if 0 in data_distrib: 
                    n = np.sum(data_distrib) + i
                    w = [data_distrib[0]/n, data_distrib[1]/n, i/n]
                else:
                    n = data_distrib[0] * i + data_distrib[1] * i + data_distrib[2] * i
                    w = [data_distrib[0] * i/n, data_distrib[1] * i/n, data_distrib[2] * i/n]
                # Features are fixed, taken from function arguments
                n_ft = np.sum(feature_distrib)
                n_ft_inf = feature_distrib[0]
                n_ft_rep = feature_distrib[1]
                # Generate dataset according to function arguments
                X, y = make_classification(n_samples=n, n_classes=len(data_distrib), n_clusters_per_class=1, 
                                           n_features=n_ft, n_informative=n_ft_inf, n_redundant=0, 
                                           n_repeated=n_ft_rep, flip_y=0, weights=w, random_state=k)
                # The following if/elif statement is here to check the distribution was correctly respected. 
                # Sometimes there is one datapoint which should be in one class but was assigned into another one.
                # This fixes this problem.
                if title == 'Datasets (a): (100, 5, i)':
                    tst = 'a'
                    if len(index)!=0 and Counter(y)[2] != index[-1] + 1:
                        if Counter(y)[0] == data_distrib[0] + 1:
                            id_y = list(y).index(0)
                            y[id_y] = 2
                    elif len(index)!=0 and Counter(y)[0] != data_distrib[0]:
                        if Counter(y)[1] == data_distrib[1] - 1:
                            id_y = list(y).index(0)
                            y[id_y] = 1
                elif title == 'Datasets (b): (50, 50, i)':
                    tst = 'b'
                    if len(index)!=0 and Counter(y)[2] != i:
                        if Counter(y)[0] == data_distrib[0] + 1:
                            id_y = list(y).index(0)
                            y[id_y] = 2
                else:
                    tst = 'c'
            else:
                # Varying features
                # Data distribution is fixed and taken from the function parameters.
                n = np.sum(data_distrib)
                w = [data_distrib[0]/n, data_distrib[1]/n, data_distrib[2]/n]
                n_ft = np.sum(feature_distrib)
                if feature_distrib[0] < 0: # (i, 0, 50-i)
                    tst = 'a'
                    n_ft_inf = i
                    n_ft_rep = 0
                    n_ft = tot_ft
                elif feature_distrib[1] < 0: # (5, i, 0)
                    if tot_ft < 0:
                        tst = 'b'
                        n_ft_inf = feature_distrib[0]
                        n_ft_rep = i
                        n_ft = n_ft_inf + i
                    else: # (50-i, i, 0)
                        tst = 'c'
                        n_ft_inf = tot_ft - i
                        n_ft_rep = i
                        n_ft = tot_ft
                # Dataset generation according to function arguments
                X, y = make_classification(n_samples=n, n_classes=len(data_distrib), n_clusters_per_class=1, 
                                           n_features=n_ft, n_informative=n_ft_inf, n_redundant=0, 
                                           n_repeated=n_ft_rep, flip_y=0, weights=w, random_state=k)
                  
            index.append(i)
            columns = ['ft' + str(j+1) for j in range(len(X[0]))]
            df = pd.DataFrame(columns=columns)
            for j in range(len(columns)):
                df[columns[j]] = X[:,j]
            df['label'] = y
            tmp = df.drop(columns=['label'])
            # metrics calculation for the dataset
            IR, adjIR, C1, C2, ID, LRID, IF, SIMBA, eval_score5cv = get_metrics_eval_score(tmp, df['label'], onlyMetrics=not corr, eval_method=eval_method)

            eval_scores.append(eval_score5cv)
            all_IR.append(IR)
            all_adjIR.append(adjIR)
            all_C1.append(C1)
            all_C2.append(C2)
            all_ID.append(ID)
            all_LRID.append(LRID)
            all_IF.append(IF)
            all_SIMBA.append(SIMBA)
            
        if len(avg_all_IR) == 0:
            avg_all_IR = all_IR
            avg_all_adjIR = all_adjIR
            avg_all_C1 = all_C1
            avg_all_C2 = all_C2
            avg_all_ID = all_ID
            avg_all_LRID = all_LRID
            avg_all_IF = all_IF
            avg_all_SIMBA = all_SIMBA
            avg_eval_scores = eval_scores
            avg_index = index
        else:
            avg_all_IR = [a + b for (a,b) in zip(all_IR, avg_all_IR)]
            avg_all_adjIR = [a + b for (a,b) in zip(all_adjIR, avg_all_adjIR)]
            avg_all_C1 = [a + b for (a,b) in zip(all_C1, avg_all_C1)]
            avg_all_C2 = [a + b for (a,b) in zip(all_C2, avg_all_C2)]
            avg_all_ID = [a + b for (a,b) in zip(all_ID, avg_all_ID)]
            avg_all_LRID = [a + b for (a,b) in zip(all_LRID, avg_all_LRID)]
            avg_all_IF = [a + b for (a,b) in zip(all_IF, avg_all_IF)]
            avg_all_SIMBA = [a + b for (a,b) in zip(all_SIMBA, avg_all_SIMBA)]
            avg_eval_scores = [a + b for (a,b) in zip(eval_scores, avg_eval_scores)]
            avg_index = [a + b for (a,b) in zip(index, avg_index)]
        eval_scores = np.array(eval_scores).reshape(-1, 1)
        all_IR = np.array(all_IR)
        all_adjIR = np.array(all_adjIR)
        all_C1 = np.array(all_C1)
        all_C2 = np.array(all_C2)
        all_ID = np.array(all_ID)
        all_LRID = np.array(all_LRID)
        all_IF = np.array(all_IF)
        all_SIMBA = np.array(all_SIMBA)
        if corr:
            # PCC Correlation with f1/gmean score
            avg_corr_IR, l_corr_IR = add_correlation(eval_scores, all_IR, avg_corr_IR, l_corr_IR)
            avg_corr_adjIR, l_corr_adjIR = add_correlation(eval_scores, all_adjIR, avg_corr_adjIR, l_corr_adjIR)
            avg_corr_C1, l_corr_C1 = add_correlation(eval_scores, all_C1, avg_corr_C1, l_corr_C1)
            avg_corr_C2, l_corr_C2 = add_correlation(eval_scores, all_C2, avg_corr_C2, l_corr_C2)
            avg_corr_ID, l_corr_ID = add_correlation(eval_scores, all_ID, avg_corr_ID, l_corr_ID)
            avg_corr_LRID, l_corr_LRID = add_correlation(eval_scores, all_LRID, avg_corr_LRID, l_corr_LRID)
            avg_corr_IF, l_corr_IF = add_correlation(eval_scores, all_IF, avg_corr_IF, l_corr_IF)
            avg_corr_SIMBA, l_corr_SIMBA = add_correlation(eval_scores, all_SIMBA, avg_corr_SIMBA, l_corr_SIMBA)
            # SRCC Correlation with f1/gmean score
            avg_corr_IR_srcc += spearmanr(eval_scores, all_IR).statistic
            avg_corr_adjIR_srcc += spearmanr(eval_scores, all_adjIR).statistic
            avg_corr_C1_srcc += spearmanr(eval_scores, all_C1).statistic
            avg_corr_C2_srcc += spearmanr(eval_scores, all_C2).statistic
            avg_corr_ID_srcc += spearmanr(eval_scores, all_ID).statistic
            avg_corr_LRID_srcc += spearmanr(eval_scores, all_LRID).statistic
            avg_corr_IF_srcc += spearmanr(eval_scores, all_IF).statistic
            avg_corr_SIMBA_srcc += spearmanr(eval_scores, all_SIMBA).statistic
            l_corr_IR_srcc.append(spearmanr(eval_scores, all_IR).statistic)
            l_corr_adjIR_srcc.append(spearmanr(eval_scores, all_adjIR).statistic)
            l_corr_C1_srcc.append(spearmanr(eval_scores, all_C1).statistic)
            l_corr_C2_srcc.append(spearmanr(eval_scores, all_C2).statistic)
            l_corr_ID_srcc.append(spearmanr(eval_scores, all_ID).statistic)
            l_corr_LRID_srcc.append(spearmanr(eval_scores, all_LRID).statistic)
            l_corr_IF_srcc.append(spearmanr(eval_scores, all_IF).statistic)
            l_corr_SIMBA_srcc.append(spearmanr(eval_scores, all_SIMBA).statistic)

    all_IR = [i/N for i in avg_all_IR]
    all_adjIR = [i/N for i in avg_all_adjIR]
    all_C1 = [i/N for i in avg_all_C1]
    all_C2 = [i/N for i in avg_all_C2]
    all_ID = [i/N for i in avg_all_ID]
    all_LRID = [i/N for i in avg_all_LRID]
    all_IF = [i/N for i in avg_all_IF]
    all_SIMBA = [i/N for i in avg_all_SIMBA]
    index = [i/N for i in avg_index]
    eval_scores = [eval_score/N for eval_score in avg_eval_scores]

    res = pd.DataFrame(columns=['index', 'IR', 'AdjIR', 'C1', 'C2', 'ID', 'LRID', 'IF', 'SIMBA', eval_method])
    res_corr = pd.DataFrame(columns=[
                        'PCC_IR', 'PCC_AdjIR', 'PCC_C1', 'PCC_C2', 'PCC_ID', 'PCC_LRID', 'PCC_IF', 'PCC_SIMBA',
                        'SRCC_IR', 'SRCC_AdjIR', 'SRCC_C1', 'SRCC_C2', 'SRCC_ID', 'SRCC_LRID', 'SRCC_IF', 'SRCC_SIMBA'])
    res['index'] = index
    res['IR'] = all_IR
    res['AdjIR'] = all_adjIR
    res['C1'] = all_C1
    res['C2'] = all_C2
    res['ID'] = all_ID
    res['LRID'] = all_LRID
    res['IF'] = all_IF
    res['SIMBA'] = all_SIMBA
    res[eval_method] = eval_scores

    metrics = ['IR', 'AdjIR', 'C1', 'C2', 'ID', 'LRID', 'IF', 'SIMBA']
    l_metrics_pcc = [l_corr_IR, l_corr_adjIR, l_corr_C1, l_corr_C2, l_corr_ID, l_corr_LRID, l_corr_IF, l_corr_SIMBA]
    l_metrics_srcc = [l_corr_IR_srcc, l_corr_adjIR_srcc, l_corr_C1_srcc, l_corr_C2_srcc, l_corr_ID_srcc, 
                      l_corr_LRID_srcc, l_corr_IF_srcc, l_corr_SIMBA_srcc]
    for k in range(len(metrics)):
        metric = metrics[k]
        try:
            res_corr['PCC_' + metric] = [i[0] for i in l_metrics_pcc[k]]
        except:
            res_corr['PCC_IR'] = l_metrics_pcc[k] 
        res_corr['SRCC_' + metric] = l_metrics_srcc[k] 
    
    # Save results for future plotting and correlation results
    res.to_csv('./table/results/plot_synth_varying_' + varying + '_' + tst + '_' + eval_method + '_C2.csv', index=False)
    res_corr = res_corr.fillna(0)
    res_corr.to_csv('./table/results/corr_synth_varying_' + varying + '_' + tst + '_' + eval_method + '_C2.csv', index=False)
    if not corr:
        # Return metrics
        return index, all_IR, all_adjIR, all_C1, all_C2, all_ID, all_LRID, all_IF, all_SIMBA
    else: # shows correlation results
        print("PCC (mean +- sd)")
        print('Correlation ' + eval_method + ' scores & IR: ' + str(round(np.mean(l_corr_IR), 2)) + " +- " + str(round(np.std(l_corr_IR), 2)))
        print('Correlation ' + eval_method + ' scores & Adj-IR: ' + str(round(np.mean(l_corr_adjIR), 2)) + " +- " + str(round(np.std(l_corr_adjIR), 2)))
        print('Correlation ' + eval_method + ' scores & C1: ' + str(round(np.mean(l_corr_C1), 2)) + " +- " + str(round(np.std(l_corr_C1), 2)))
        print('Correlation ' + eval_method + ' scores & C2: ' + str(round(np.mean(l_corr_C2), 2)) + " +- " + str(round(np.std(l_corr_C2), 2)))
        print('Correlation ' + eval_method + ' scores & ID: ' + str(round(np.mean(l_corr_ID), 2)) + " +- " + str(round(np.std(l_corr_ID), 2)))
        print('Correlation ' + eval_method + ' scores & LRID: ' + str(round(np.mean(l_corr_LRID), 2)) + " +- " + str(round(np.std(l_corr_LRID), 2)))
        print('Correlation ' + eval_method + ' scores & IF: ' + str(round(np.mean(l_corr_IF), 2)) + " +- " + str(round(np.std(l_corr_IF), 2)))
        print('Correlation ' + eval_method + ' scores & SIMBA: ' + str(round(np.mean(l_corr_SIMBA), 2)) + " +- " + str(round(np.std(l_corr_SIMBA), 2)))
        print()
        print("SRCC (mean +- sd)")
        print('Correlation ' + eval_method + ' scores & IR: ' + str(round(np.mean(l_corr_IR_srcc), 2)) + " +- " + str(round(np.std(l_corr_IR_srcc), 2)))
        print('Correlation ' + eval_method + ' scores & Adj-IR: ' + str(round(np.mean(l_corr_adjIR_srcc), 2)) + " +- " + str(round(np.std(l_corr_adjIR_srcc), 2)))
        print('Correlation ' + eval_method + ' scores & C1: ' + str(round(np.mean(l_corr_C1_srcc), 2)) + " +- " + str(round(np.std(l_corr_C1_srcc), 2)))
        print('Correlation ' + eval_method + ' scores & C2: ' + str(round(np.mean(l_corr_C2_srcc), 2)) + " +- " + str(round(np.std(l_corr_C2_srcc), 2)))
        print('Correlation ' + eval_method + ' scores & ID: ' + str(round(np.mean(l_corr_ID_srcc), 2)) + " +- " + str(round(np.std(l_corr_ID_srcc), 2)))
        print('Correlation ' + eval_method + ' scores & LRID: ' + str(round(np.mean(l_corr_LRID_srcc), 2)) + " +- " + str(round(np.std(l_corr_LRID_srcc), 2)))
        print('Correlation ' + eval_method + ' scores & IF: ' + str(round(np.mean(l_corr_IF_srcc), 2)) + " +- " + str(round(np.std(l_corr_IF_srcc), 2)))
        print('Correlation ' + eval_method + ' scores & SIMBA: ' + str(round(np.mean(l_corr_SIMBA_srcc), 2)) + " +- " + str(round(np.std(l_corr_SIMBA_srcc), 2)))

def run_synthdata_plot(N, varying="data", eval_method="f1", from_file=False): 
    """
    Creates and runs all the tests on synthetic datasets, but only on the metrics (no correlations). This function
    will generate plots for each imbalance metrics and show their trend on the three data distribution variations
    (if varying=="data") or the three features distribution variations (if varying!="data").
    @N: integer, number of random runs. For each run, a different random state is picked. 
    @varying: String, defines what is varying between "data" or features. If == "data", data distribution is varying.
              If != "data", it is considered that the features are varying.
    @eval_method: String, defines the classification evaluation metric: either "f1" for f1-score or gmean
    @from_file: Boolean, default=False. When False, the synthetic data analysis is run with run_synthdata(). When True,
                the run_synthdata() function has already been run and the plots are generated from the files saved with
                this function.
    """
    if varying == "data": # varying data distribution
        x_lim=[0, 105]
        if not from_file:
            index_a, all_IR_a, all_adjIR_a, all_C1_a, all_C2_a, all_ID_a, all_LRID_a, all_IF_Collision_a, all_new_metric_a =\
            run_synthdata(N = N, data_distrib = [100, 5, 0], feature_distrib = [2, 0, 8], tot_ft = 10, 
                        i_range=[5, 101, 1], title='Datasets (a): (100, 5, i)', corr=False, varying="data", 
                        eval_method=eval_method)
            index_b, all_IR_b, all_adjIR_b, all_C1_b, all_C2_b, all_ID_b, all_LRID_b, all_IF_Collision_b, all_new_metric_b =\
            run_synthdata(N = N, data_distrib = [50, 50, 0], feature_distrib = [2, 0, 8], tot_ft = 10, 
                        i_range=[5, 101, 1], title='Datasets (b): (50, 50, i)', corr=False, varying="data",
                        eval_method=eval_method)
            index_c, all_IR_c, all_adjIR_c, all_C1_c, all_C2_c, all_ID_c, all_LRID_c, all_IF_Collision_c, all_new_metric_c =\
            run_synthdata(N = N, data_distrib = [10, 50, 100], feature_distrib = [2, 0, 8], tot_ft = 10, 
                        i_range=[5, 101, 1], title='Datasets (c): (10i, 50i, 100i)', corr=False, varying="data",
                        eval_method=eval_method)   
        else:
            file_a = pd.read_csv('./table/results/plot_synth_varying_data_a_' + eval_method + '.csv')
            index_a, all_IR_a, all_adjIR_a, all_C1_a, all_C2_a, all_ID_a, all_LRID_a, all_IF_Collision_a, all_new_metric_a =\
            get_values_plot_synth(file_a)
            file_b = pd.read_csv('./table/results/plot_synth_varying_data_b_' + eval_method + '.csv')
            index_b, all_IR_b, all_adjIR_b, all_C1_b, all_C2_b, all_ID_b, all_LRID_b, all_IF_Collision_b, all_new_metric_b =\
            get_values_plot_synth(file_b)
            file_c = pd.read_csv('./table/results/plot_synth_varying_data_c_' + eval_method + '.csv')
            index_c, all_IR_c, all_adjIR_c, all_C1_c, all_C2_c, all_ID_c, all_LRID_c, all_IF_Collision_c, all_new_metric_c =\
            get_values_plot_synth(file_c)

        metrics = ["IR", "Adj-IR", "C1", "C2", "ID", "LRID", "IF", "SIMBA"]
        # These are chosen for clearer readibility of the plots.
        lims = [[[15, 21], [0, 11], [9, 15]], # IR
                [[15, 21], [0, 11], [8.5, 15]], # Adj-IR
                [[0.3, 0.8], [0.7, 1.025], [0.7, 0.8]], # C1
                [[0.2, 0.9], [-0.01, 0.21], [0., 1.]], # C2 TODO change lims
                [[0.6, 1.8], [-1.1, 1.5], [1, 2]], # ID
                [[80, 130], [-1, 45], [0, 8000]], # LRID
                [[0.1, 0.7], [0.7, 1.025], [0.6, 0.7]], # IF
                [[4, 11], [-0.1,3.2], [1, 10]]] # SIMBA
    else: # varying feature distribution
        x_lim = [-1, 51]
        if not from_file:
            index_a, all_IR_a, all_adjIR_a, all_C1_a, all_C2_a, all_ID_a, all_LRID_a, all_IF_Collision_a, all_new_metric_a =\
            run_synthdata(N = N, data_distrib = [400, 75, 25], feature_distrib = [-1, 0, 0], tot_ft = 50, 
                        i_range=[2, 51, 1], title='Datasets with i out of 50 informative features, dataset (400, 75, 25)', 
                        corr=False, varying="ft", eval_method=eval_method)
            index_b, all_IR_b, all_adjIR_b, all_C1_b, all_C2_b, all_ID_b, all_LRID_b, all_IF_Collision_b, all_new_metric_b =\
            run_synthdata(N = N, data_distrib = [400, 75, 25], feature_distrib = [5, -1, 0], tot_ft = -1, 
                        i_range=[0, 46, 1], title='Datasets with 5+i features, 5 informative, and i redundant, dataset (400, 75, 25)', 
                        corr=False, varying="ft", eval_method=eval_method)
            index_c, all_IR_c, all_adjIR_c, all_C1_c, all_C2_c, all_ID_c, all_LRID_c, all_IF_Collision_c, all_new_metric_c =\
            run_synthdata(N = N, data_distrib = [400, 75, 25], feature_distrib = [0, -1, 0], tot_ft = 50, 
                        i_range=[0, 49, 1], title='Datasets with 50 features, 50-i informative, and i redundant, dataset (400, 75, 25)', 
                        corr=False, varying="ft", eval_method=eval_method)
        else:
            file_a = pd.read_csv('./table/results/plot_synth_varying_ft_a_' + eval_method + '.csv')
            index_a, all_IR_a, all_adjIR_a, all_C1_a, all_C2_a, all_ID_a, all_LRID_a, all_IF_Collision_a, all_new_metric_a =\
            get_values_plot_synth(file_a)
            file_b = pd.read_csv('./table/results/plot_synth_varying_ft_b_' + eval_method + '.csv')
            index_b, all_IR_b, all_adjIR_b, all_C1_b, all_C2_b, all_ID_b, all_LRID_b, all_IF_Collision_b, all_new_metric_b =\
            get_values_plot_synth(file_b)
            file_c = pd.read_csv('./table/results/plot_synth_varying_ft_c_' + eval_method + '.csv')
            index_c, all_IR_c, all_adjIR_c, all_C1_c, all_C2_c, all_ID_c, all_LRID_c, all_IF_Collision_c, all_new_metric_c =\
            get_values_plot_synth(file_c)
        metrics = ["IR", "Adj-IR", "C1", "C2", "ID", "LRID", "IF", "SIMBA"]
        # Only Adj-IR and SIMBA vary when data distribution is fixed and features distribution varies. Indeed, all
        # the other metrics remain constant as they do not take features into account at all. 
        lims = [[],
                [[11, 14.5], [10.5, 14.5], [10.3, 11.6]], # Adj-IR
                [],
                [],
                [],
                [],
                [],
                [[10, 18], [3, 8], [1, 11]]] # SIMBA
        
    all_metrics_a = [all_IR_a, all_adjIR_a, all_C1_a, all_C2_a, all_ID_a, all_LRID_a, all_IF_Collision_a, all_new_metric_a]
    all_metrics_b = [all_IR_b, all_adjIR_b, all_C1_b, all_C2_b, all_ID_b, all_LRID_b, all_IF_Collision_b, all_new_metric_b]
    all_metrics_c = [all_IR_c, all_adjIR_c, all_C1_c, all_C2_c, all_ID_c, all_LRID_c, all_IF_Collision_c, all_new_metric_c]

    # Plots generation
    for i in range(len(metrics)): # one plot per imbalance metric
        metric = metrics[i]
        if varying == "data" or metric == "Adj-IR" or metric == "SIMBA":
            lim = lims[i]
            all_metric_a = all_metrics_a[i]
            all_metric_b = all_metrics_b[i]
            all_metric_c = all_metrics_c[i]
            _, ax_a = plt.subplots(figsize=(5,4))
            ax_b = ax_a.twinx()
            ax_c = ax_a.twinx()
            if varying == 'data':
                ax_a.set_xticks([5,50,100])
                ax_a.set_xlim(x_lim[0], x_lim[1])
                if min(all_metric_a) == max(all_metric_a):
                    ax_a.set_yticks([min(all_metric_a)])
                else:
                    if metric == "Adj-IR":
                        ax_a.set_yticks([np.mean(all_metric_a)])
                    else:
                        step = (max(all_metric_a) - min(all_metric_a))/5
                        ax_a.set_yticks([min(all_metric_a) + i*step for i in range(0,5)] + [max(all_metric_a)])
                if min(all_metric_b) == max(all_metric_b):
                    ax_b.set_yticks([min(all_metric_b)])
                else:
                    step = (max(all_metric_b) - min(all_metric_b))/5
                    ax_b.set_yticks([min(all_metric_b) + i*step for i in range(0,5)] + [max(all_metric_b)])
                if min(all_metric_c) == max(all_metric_c):
                        ax_c.set_yticks([min(all_metric_c)])
                else:
                    if metric == "Adj-IR":
                        ax_c.set_yticks([np.mean(all_metric_c)])
                    elif metric == "SIMBA":
                        ax_c.set_yticks([min(all_metric_c), max(all_metric_c)])
                    else:
                        step = (max(all_metric_c) - min(all_metric_c))/5
                        ax_c.set_yticks([min(all_metric_c) + i*step for i in range(0,5)] + [max(all_metric_c)])
            else:
                ax_a.set_xticks([0, 5, 45, 50])
                ax_a.set_xlim(x_lim[0], x_lim[1])
                if min(all_metric_a) == max(all_metric_a):
                    ax_a.set_yticks([min(all_metric_a)])
                else:
                    step = (max(all_metric_a) - min(all_metric_a))/5
                    ax_a.set_yticks([min(all_metric_a) + i*step for i in range(0,5)] + [max(all_metric_a)])
                if min(all_metric_b) == max(all_metric_b):
                    ax_b.set_yticks([min(all_metric_b)])
                else:
                    if metric == "SIMBA":
                        ax_b.set_yticks([np.mean(all_metric_b)])
                    else:
                        step = (max(all_metric_b) - min(all_metric_b))/5
                        ax_b.set_yticks([min(all_metric_b) + i*step for i in range(0,5)] + [max(all_metric_b)])
                if min(all_metric_c) == max(all_metric_c):
                        ax_c.set_yticks([min(all_metric_c)])
                else:
                    step = (max(all_metric_c) - min(all_metric_c))/5
                    ax_c.set_yticks([min(all_metric_c) + i*step for i in range(0,5)] + [max(all_metric_c)])
            if metric == 'IF' or metric == 'C1' or metric == 'C2':
                ax_a.set_yticklabels([round(a, 2) for a in ax_a.get_yticks()])
                ax_b.set_yticklabels([round(a, 2) for a in ax_b.get_yticks()])
                ax_c.set_yticklabels([round(a, 2) for a in ax_c.get_yticks()])
            else:
                ax_a.set_yticklabels([round(a, 1) for a in ax_a.get_yticks()])
                ax_b.set_yticklabels([round(a, 1) for a in ax_b.get_yticks()])
                ax_c.set_yticklabels([round(a, 1) for a in ax_c.get_yticks()])
            ax_a.set_ylim(lim[0][0], lim[0][1])
            ax_b.set_ylim(lim[1][0], lim[1][1])
            ax_c.set_ylim(lim[2][0], lim[2][1])
            ax_a.set_xlabel("i", fontsize=17.5)
            ax_a.set_ylabel(metric, fontsize=17.5)
            if varying == "data":
                color1 = '#882255'
                color2 = '#44AA99'
                color3 = '#88CCEE'
                p1 = ax_a.plot(index_a, all_metric_a, color=color1, label="(100, 5, i)", linewidth=7)
                p2 = ax_b.plot(index_b, all_metric_b, color=color2, label="(50, 50, i)", linewidth=7)
                p3 = ax_c.plot(index_c, all_metric_c, color=color3, label="(10i, 50i, 100i)", linewidth=7)
            else:
                color1 = '#332288'
                color2 = '#117733'
                color3 = '#CC6677'               
                p1 = ax_a.plot(index_a, all_metric_a, color=color1, label="(i, 0, 50-i)", linewidth=7)
                p2 = ax_b.plot(index_b, all_metric_b, color=color2, label="(5, i, 0)", linewidth=7)
                p3 = ax_c.plot(index_c, all_metric_c, color=color3, label="(50-i, i, 0)", linewidth=7)
            
            ax_a.tick_params(axis='x', labelsize=14.5)
            ax_a.tick_params(axis='y', colors=p1[0].get_color(), labelsize=14.5)
            ax_b.tick_params(axis='y', colors=p2[0].get_color(), labelsize=14.5)
            ax_c.tick_params(axis='y', colors=p3[0].get_color(), labelsize=14.5)
            ax_c.spines['right'].set_position(('outward', 40))
            plt.tight_layout()
            plt.show()

def run_realdata(classif="SVM", eval_method="f1"):
    """
    Runs the tests on all real datasets (contained in data/ folder). Metrics and f1-scores/gmeans are computed for each
    dataset. 
    @eval_method: String, defines the classification evaluation metric: either "f1" for f1-score or gmean
    """

    # list of all datasets: multiclass datasets are in the l_datasets and binary ones are in l_datasets_bin
    # in order for us to be able to separate them for the plots.
    l_datasets = ['abalone', 'balance', 'cardio10', 'cardio3', 'chess', 
                'connect-4', 'contraceptive', 'dermatology', 'drybean', 'ecoli', 
                'glass', 'hayes-roth', 'knowledge', 'landsat', 'lenses', 
                'loc_build', 'loc_floor', 'lymphography', 'new-thyroid', 'obesity', 
                'pageblocks', 'penbased', 'room', 'shuttle', 'soybean', 
                'steel', 'student', 'theorem', 'thyroid', 'wallfollowing', 
                'webphishing', 'wholesale', 'wine', 'wine-quality', 'yeast']
    l_datasets_bin = ['abalone-20_vs_8-9-10', 'abalone9-18', 'adult', 'banknote', 'bankruptcy',
                    'breastcancer', 'cleveland-0_vs_4', 'credit', 'dermatology-6', 'ecoli-0-1-4-7_vs_2-3-5-6', 
                    'ecoli1', 'glass-0-1-5_vs_2', 'glass2', 'glass4', 'glass5', 
                    'glass6', 'htru2', 'ionosphere', 'led7digit-0-2-4-5-6-7-8-9_vs_1', 'new-thyroid1', 
                    'page-blocks0', 'poker-9_vs_7', 'purchase', 'segment0', 'skin',
                    'spambase', 'spect-heart', 'vehicle0', 'vowel0', 'winequality-red-3_vs_5', 
                    'wisconsin', 'yeast-1-2-8-9_vs_7', 'yeast3', 'yeast4', 'yeast5']
    l_datasets = l_datasets + l_datasets_bin
    # all datasets mentioned in the lists above have 'label' and 'id' respectively as ID and labels columns' names
    l_col_label = ["label"] * len(l_datasets)
    l_col_ID = [["id"]] * len(l_datasets)
   
    all_IR = []
    all_adjIR = []
    all_C1 = []
    all_C2 = []
    all_ID = []
    all_LRID = []
    all_IF = []
    all_SIMBA = []
    eval_scores = []
    binary = []
    columns = ['DS', 'binary', 'IR', 'AdjIR', 'C1', 'C2', 'ID', 'LRID', 'IF', 'SIMBA', eval_method]
    res = pd.DataFrame(columns=columns)
    for i in range(len(l_datasets)):
        dataset = l_datasets[i]
        if dataset in l_datasets_bin:
            binary.append(1) # True
        else:
            binary.append(0) # False, multiclass
        print("DATASET: " + str(dataset))
        col_label = l_col_label[i]
        col_ID = l_col_ID[i]

        path = "./data/" + dataset + "/"
        X = pd.read_csv(path + "df_allfeatures.csv") # features data
        y = pd.read_csv(path + "labels.csv") # labels data

        df = X.drop(columns=col_ID) # drops id column
        df[col_label] = y[col_label] # adds label data in df
        tmp = df.drop(columns=[col_label]) # drops the labels

        # metrics and f1-scores/gmean calculation
        IR, adjIR, C1, C2, ID, LRID, IF, SIMBA, eval_scores5cv = get_metrics_eval_score(tmp, df[col_label], 
                                                                                        classif=classif, 
                                                                                        onlyMetrics=False, 
                                                                                        eval_method=eval_method)
        all_IR.append(IR)
        all_adjIR.append(adjIR)
        all_C1.append(C1)
        all_C2.append(C2)
        all_ID.append(ID)
        all_LRID.append(LRID)
        all_IF.append(IF)
        all_SIMBA.append(SIMBA)
        eval_scores.append(eval_scores5cv)
        
    res['DS'] = l_datasets
    res['binary'] = binary
    res['IR'] = all_IR
    res['AdjIR'] = all_adjIR
    res['C1'] = all_C1
    res['C2'] = all_C2
    res['ID'] = all_ID
    res['LRID'] = all_LRID
    res['IF'] = all_IF
    res['SIMBA'] = all_SIMBA
    res[eval_method] = eval_scores
    
    res.to_csv('./table/results/real_' + str(classif) + '_' + str(eval_method) + 'C2.csv', index=False)
    print('Results saved!')
    eval_scores = np.array(eval_scores).reshape(-1, 1)
    
    print("PCC")
    print('Correlation ' + eval_method + ' scores & IR: ' + str(r_regression(eval_scores, all_IR)))
    print('Correlation ' + eval_method + ' scores & adjIR: ' + str(r_regression(eval_scores, all_adjIR)))
    print('Correlation ' + eval_method + ' scores & C1: ' + str(r_regression(eval_scores, all_C1)))
    print('Correlation ' + eval_method + ' scores & C2: ' + str(r_regression(eval_scores, all_C2)))
    print('Correlation ' + eval_method + ' scores & ID: ' + str(r_regression(eval_scores, all_ID)))
    print('Correlation ' + eval_method + ' scores & LRID: ' + str(r_regression(eval_scores, all_LRID)))
    print('Correlation ' + eval_method + ' scores & IF: ' + str(r_regression(eval_scores, all_IF)))
    print('Correlation ' + eval_method + ' scores & SIMBA: ' + str(r_regression(eval_scores, all_SIMBA)))

    print("SRCC")
    print('Correlation ' + eval_method + ' scores & IR: ' + str(spearmanr(eval_scores, all_IR)))
    print('Correlation ' + eval_method + ' scores & adjIR: ' + str(spearmanr(eval_scores, all_adjIR)))
    print('Correlation ' + eval_method + ' scores & C1: ' + str(spearmanr(eval_scores, all_C1)))
    print('Correlation ' + eval_method + ' scores & C2: ' + str(spearmanr(eval_scores, all_C2)))
    print('Correlation ' + eval_method + ' scores & ID: ' + str(spearmanr(eval_scores, all_ID)))
    print('Correlation ' + eval_method + ' scores & LRID: ' + str(spearmanr(eval_scores, all_LRID)))
    print('Correlation ' + eval_method + ' scores & IF: ' + str(spearmanr(eval_scores, all_IF)))
    print('Correlation ' + eval_method + ' scores & SIMBA: ' + str(spearmanr(eval_scores, all_SIMBA)))

def run_real_ablation(classif="SVM", eval_method="f1"):
    """
    Runs the ablation study on all real datasets (contained in data/ folder). 4 versions of SIMBA (whole, without 
    normalization, without feature importance, and without feature redundancy) and f1-scores/gmeans are computed 
    for each dataset. PCC & SRCC between SIMBA scores (in all versions) and f1-score/gmean for all datasets are 
    generated.
    @classif: (optional) String, default="SVM", chosen classifier, can take the following values: "SVM", "LDA", 
              "kNN", "RF", and "MLP". 
              [WARNING!] any other value will lead to an error and the program will be stopped.
    @eval_method: String, defines the classification evaluation metric: either "f1" for f1-score or gmean
    """
    l_datasets = ['abalone', 'balance', 'cardio10', 'cardio3', 'chess', 
                'connect-4', 'contraceptive', 'dermatology', 'drybean', 'ecoli', 
                'glass', 'hayes-roth', 'knowledge', 'landsat', 'lenses', 
                'loc_build', 'loc_floor', 'lymphography', 'new-thyroid', 'obesity', 
                'pageblocks', 'penbased', 'room', 'shuttle', 'soybean', 
                'steel', 'student', 'theorem', 'thyroid', 'wallfollowing', 
                'webphishing', 'wholesale', 'wine', 'wine-quality', 'yeast']
    l_datasets_bin = ['abalone-20_vs_8-9-10', 'abalone9-18', 'adult', 'banknote', 'bankruptcy',
                    'breastcancer', 'cleveland-0_vs_4', 'credit', 'dermatology-6', 'ecoli-0-1-4-7_vs_2-3-5-6', 
                    'ecoli1', 'glass-0-1-5_vs_2', 'glass2', 'glass4', 'glass5', 
                    'glass6', 'htru2', 'ionosphere', 'led7digit-0-2-4-5-6-7-8-9_vs_1', 'new-thyroid1', 
                    'page-blocks0', 'poker-9_vs_7', 'purchase', 'segment0', 'skin',
                    'spambase', 'spect-heart', 'vehicle0', 'vowel0', 'winequality-red-3_vs_5', 
                    'wisconsin', 'yeast-1-2-8-9_vs_7', 'yeast3', 'yeast4', 'yeast5']
    l_datasets = l_datasets + l_datasets_bin
    # all datasets mentioned in the lists above have 'label' and 'id' respectively as ID and labels columns' names
    l_col_label = ["label"] * len(l_datasets)
    l_col_ID = [["id"]] * len(l_datasets)
   
    all_SIMBA_whole = []
    all_SIMBA_nonorm = []
    all_SIMBA_noftimp = []
    all_SIMBA_noftred = []
    eval_scores = []

    res_ablation = pd.DataFrame(columns=['DS', 'SIMBA_whole', 'SIMBA_nonorm', 'SIMBA_noftimp', 'SIMBA_noftred', eval_method])
    for i in range(len(l_datasets)):
        
        dataset = l_datasets[i]
        print("DATASET: " + str(dataset))
        col_label = l_col_label[i]
        col_ID = l_col_ID[i]

        path = "./data/" + dataset + "/"
        X = pd.read_csv(path + "df_allfeatures.csv") # features data
        y = pd.read_csv(path + "labels.csv") # labels data

        df = X.drop(columns=col_ID) # drops id column
        df[col_label] = y[col_label] # adds label data in df
        tmp = df.drop(columns=[col_label]) # drops the labels

        # metrics and f1-scores/gmean calculation
        SIMBA_whole, SIMBA_nonorm, SIMBA_noftimp, SIMBA_noftred, eval_score = get_SIMBA_ablation_eval_score(tmp, df[col_label], classif=classif)
        all_SIMBA_whole.append(SIMBA_whole)
        all_SIMBA_nonorm.append(SIMBA_nonorm)
        all_SIMBA_noftimp.append(SIMBA_noftimp)
        all_SIMBA_noftred.append(SIMBA_noftred)
        eval_scores.append(eval_score)
        
    res_ablation['DS'] = l_datasets
    res_ablation['SIMBA_whole'] = all_SIMBA_whole
    res_ablation['SIMBA_nonorm'] = all_SIMBA_nonorm
    res_ablation['SIMBA_noftimp'] = all_SIMBA_noftimp
    res_ablation['SIMBA_noftred'] = all_SIMBA_noftred
    res_ablation[eval_method] = eval_scores
    res_ablation.to_csv('./table/results/ablation_real_' + str(classif) + '_' + str(eval) + '.csv', index=False)
    print('Results saved!')
    eval_scores = np.array(eval_scores).reshape(-1, 1)
    
    print("PCC")
    print('Correlation ' + eval_method + ' scores & SIMBA (whole): ' + str(r_regression(eval_scores, all_SIMBA_whole)))
    print('Correlation ' + eval_method + ' scores & SIMBA (nonorm): ' + str(r_regression(eval_scores, all_SIMBA_nonorm)))
    print('Correlation ' + eval_method + ' scores & SIMBA (noftimp): ' + str(r_regression(eval_scores, all_SIMBA_noftimp)))
    print('Correlation ' + eval_method + ' scores & SIMBA (noftred): ' + str(r_regression(eval_scores, all_SIMBA_noftred)))

    print("SRCC")
    print('Correlation ' + eval_method + ' scores & SIMBA (whole): ' + str(spearmanr(eval_scores, all_SIMBA_whole)))
    print('Correlation ' + eval_method + ' scores & SIMBA (nonorm): ' + str(spearmanr(eval_scores, all_SIMBA_nonorm)))
    print('Correlation ' + eval_method + ' scores & SIMBA (noftimp): ' + str(spearmanr(eval_scores, all_SIMBA_noftimp)))
    print('Correlation ' + eval_method + ' scores & SIMBA (noftred): ' + str(spearmanr(eval_scores, all_SIMBA_noftred)))

def plot_real_ablation(classif="SVM", eval_method="f1"):
    """
    Plots results of the ablation study that was run with the run_real_ablation function.
    [WARNING] You must run the run_real_ablation function with the same values for classif and eval_method first!
    @classif: (optional) String, default="SVM", chosen classifier, can take the following values: "SVM", "LDA", 
              "kNN", "RF", and "MLP". 
              [WARNING!] any other value will lead to an error and the program will be stopped.
    @eval_method: String, defines the classification evaluation metric: either "f1" for f1-score or gmean
    """
    file = './table/results/ablation_real_' + classif + '_' + eval_method + '.csv'
    results = pd.read_csv(file)
    res_whole = results['SIMBA_whole'].to_list()
    res_nonorm = results['SIMBA_nonorm'].to_list()
    res_noftimp = results['SIMBA_noftimp'].to_list()
    res_noftred = results['SIMBA_noftred'].to_list()
    classif_performance = results['f1'].to_list()
    classif_performance = np.array(classif_performance).reshape(-1, 1)
    # correlation calculation for all datasets, only binary datasets, only multiclass datasets
    corr_whole = np.round(r_regression(classif_performance, res_whole), 4)
    corr_nonorm = np.round(r_regression(classif_performance, res_nonorm), 4)
    corr_noftimp = np.round(r_regression(classif_performance, res_noftimp), 4)
    corr_noftred = np.round(r_regression(classif_performance, res_noftred), 4)
    corr_whole_srcc = np.round(spearmanr(classif_performance, res_whole)[0], 4)
    corr_nonorm_srcc = np.round(spearmanr(classif_performance, res_nonorm)[0], 4)
    corr_noftimp_srcc = np.round(spearmanr(classif_performance, res_noftimp)[0], 4)
    corr_noftred_srcc = np.round(spearmanr(classif_performance, res_noftred)[0], 4)
    #  plot all (show ablation study, 2x2 format)
    sns.set(style="darkgrid")
    fig, axs = plt.subplots(2, 2, figsize=(7.2, 6))
    plot_subplot(axs[0, 0], res_whole, classif_performance, 'SIMBA', corr_whole[0], corr_whole_srcc)
    plot_subplot(axs[0, 1], res_nonorm, classif_performance, 'SIMBA without normalization', corr_nonorm[0], corr_nonorm_srcc)
    plot_subplot(axs[1, 0], res_noftimp, classif_performance, 'SIMBA without feature importance', corr_noftimp[0], corr_noftimp_srcc)
    plot_subplot(axs[1, 1], res_noftred, classif_performance, 'SIMBA without feature redundancy', corr_noftred[0], corr_noftred_srcc)

    x_label = 'SIMBA score'
    y_label = 'f1-score'

    for i, ax in enumerate(axs.flat):
        ax.grid(True)
        if i in [2, 3]:  # bottom row
            ax.set_xlabel(x_label)
        if i in [0, 2]:  # left column
            ax.set_ylabel(y_label)
    plt.subplots_adjust(hspace=0.3)
    plt.subplots_adjust(wspace=0.3)
    plt.tight_layout()  
    plt.show()

def plot_subplot(ax, x, y, title, pcc, srcc):
    # function used in plot_real_ablation to plot each subplot with regression line!
    sns.regplot(x=x, y=y, ax=ax, scatter_kws={"s": 50}) 
    ax.set_title(title)
    ax.set_ylim(0, 1)  # y-axis range fixed from 0 to 1 as it is the f1-score or g-mean
    # non-overlapping annotations
    x_pos = min(x)
    ax.annotate(f'PCC = {pcc:.2f}', xy=(x_pos, 0.12), fontsize=12, weight='bold')
    ax.annotate(f'SRCC = {srcc:.2f}', xy=(x_pos, 0.05), fontsize=12, weight='bold')

def get_values_plot_synth(file_data):
    # Function that reads saved file with results, used in the run_synthdata_plot function
    index = file_data['index']
    all_IR = file_data['IR'] 
    all_adjIR = file_data['AdjIR']
    all_C1 = file_data['C1']
    all_C2 = file_data['C2'] 
    all_ID = file_data['ID'] 
    all_LRID = file_data['LRID']
    all_IF_Collision = file_data['IF']
    all_new_metric = file_data['SIMBA']
    return index, all_IR, all_adjIR, all_C1, all_C2, all_ID, all_LRID, all_IF_Collision, all_new_metric 
    
def show_realdata_res(classif, eval_method):
    """
    Function that displays the PCC and SRCC results on the real datasets and plots the regression lines.
    Plots of the PCC between metric and f1-score/gmean for all datasets, only multiclass datasets, and only
    binary datasets are generated for each imbalance metric (IR, adjIR, C1, C2, ID, LRID, IF, and SIMBA).
    [WARNING] You must run the run_realdata function with the same values for classif and eval_method first!
    @classif: (optional) String, default="SVM", chosen classifier, can take the following values: "SVM", "LDA", 
              "kNN", "RF", and "MLP". 
              [WARNING!] any other value will lead to an error and the program will be stopped.
    @eval_method: String, defines the classification evaluation metric used : either "f1" for f1-score or gmean.
    """
    print('Classifier: ' + classif)
    file = './table/results/real_' + str(classif) + '_' + str(eval_method) + '.csv'
    results = pd.read_csv(file)
    metrics = ["IR", "AdjIR", "C1", "C2", "ID", "LRID", "IF", "SIMBA"]
    eval_scores = results[eval_method].to_list()
    eval_scores_binary = results[results['binary'] == 1][eval_method].to_list()
    eval_scores_multiclass = results[results['binary'] == 0][eval_method].to_list()
    eval_scores = np.array(eval_scores).reshape(-1, 1)
    eval_scores_binary = np.array(eval_scores_binary).reshape(-1, 1)
    eval_scores_multiclass = np.array(eval_scores_multiclass).reshape(-1, 1)
    for metric in metrics:
        all_values = results[metric].to_list()
        all_values_binary = results[results['binary'] == 1][metric].to_list()
        all_values_multiclass = results[results['binary'] == 0][metric].to_list()
        tmp = pd.DataFrame(columns=[eval_method, metric, "DS"])
        tmp[eval_method] = [eval_scores[i][0] for i in range(len(eval_scores))] + [eval_scores_binary[i][0] for i in range(len(eval_scores_binary))] + [eval_scores_multiclass[i][0] for i in range(len(eval_scores_multiclass))]
        # correlation calculation for all datasets, only binary datasets, only multiclass datasets
        corr_all = np.round(r_regression(eval_scores, all_values), 4)[0]
        p_val_all = np.round(f_regression(eval_scores, np.array(all_values))[1][0],5)
        corr_bin = np.round(r_regression(eval_scores_binary, all_values_binary), 4)[0]
        p_val_bin = np.round(f_regression(eval_scores_binary, np.array(all_values_binary))[1][0],5)
        corr_multi = np.round(r_regression(eval_scores_multiclass, all_values_multiclass), 4)[0]
        p_val_multi = np.round(f_regression(eval_scores_multiclass, np.array(all_values_multiclass))[1][0],5)
        corr_all_srcc = np.round(spearmanr(eval_scores, all_values)[0], 2)
        p_val_all_srcc = np.round(spearmanr(eval_scores, all_values)[1], 5)
        corr_bin_srcc = np.round(spearmanr(eval_scores_binary, all_values_binary)[0], 4)
        p_val_bin_srcc = np.round(spearmanr(eval_scores_binary, all_values_binary)[1], 5)
        corr_multi_srcc = np.round(spearmanr(eval_scores_multiclass, all_values_multiclass)[0], 4)
        p_val_multi_srcc = np.round(spearmanr(eval_scores_multiclass, all_values_multiclass)[1], 5)
        print('PCC')
        print('Correlation ' + eval_method + ' scores & ' + metric + ': ' + str(corr_all) + ', p-value: ' + str(p_val_all))
        print('BINARY Correlation ' + eval_method + ' scores & ' + metric + ': ' + str(corr_bin) + ', p-value: ' + str(p_val_bin))
        print('MULTICLASS Correlation ' + eval_method + ' scores & ' + metric + ': ' + str(corr_multi) + ', p-value: ' + str(p_val_multi))
        print("SRCC")
        print('Correlation ' + eval_method + ' scores & ' + metric + ': ' + str(corr_all_srcc) + ', p-value: ' + str(p_val_all_srcc) )
        print('BINARY Correlation ' + eval_method + ' scores & ' + metric + ': ' + str(corr_bin_srcc) + ', p-value: ' + str(p_val_bin_srcc))
        print('MULTICLASS Correlation ' + eval_method + ' scores & ' + metric + ': ' + str(corr_multi_srcc) + ', p-value: ' + str(p_val_multi_srcc))
        input()
        tmp[metric] = all_values + all_values_binary + all_values_multiclass
        tmp["DS"] = ["All"] * len(all_values) + ["Binary"] * len(all_values_binary) + ["Multiclass"] * len(all_values_multiclass)
        #  Plot all
        sns.set(font_scale=1.5)
        g = sns.lmplot(
            data=tmp, x=metric, y=eval_method,
            col="DS", col_wrap=3,
            facet_kws=dict(sharex=False, sharey=True),
            height=5, aspect=0.85, scatter_kws={"s": 80}
        )
        g.set(ylim=(0, 1.05))
        g.set_titles('{col_name}', size=15)
        ax = g.facet_axis(0, 0)
        ax.annotate("PCC="+str(corr_all), [min(all_values), 0.075], fontsize=13, weight='bold')
        ax.annotate("SRCC="+str(corr_all_srcc), [min(all_values), 0.025], fontsize=13, weight='bold')
        ax.xaxis.get_label().set_fontsize(13)
        ax.yaxis.get_label().set_fontsize(13)
        ax.tick_params(axis='x', labelsize=12)
        ax.tick_params(axis='y', labelsize=12)
        ax = g.facet_axis(0, 1)
        ax.annotate("PCC="+str(corr_bin), [min(all_values_binary), 0.075], fontsize=13, weight='bold')
        ax.annotate("SRCC="+str(corr_bin_srcc), [min(all_values_binary), 0.025], fontsize=13, weight='bold')
        ax.xaxis.get_label().set_fontsize(13)
        ax.yaxis.get_label().set_fontsize(13)
        ax.tick_params(axis='x', labelsize=12)
        ax.tick_params(axis='y', labelsize=12)
        ax = g.facet_axis(0, 2)
        ax.annotate("PCC="+str(corr_multi), [min(all_values_multiclass), 0.075], fontsize=13, weight='bold')
        ax.annotate("SRCC="+str(corr_multi_srcc), [min(all_values_multiclass), 0.025], fontsize=13, weight='bold')
        ax.xaxis.get_label().set_fontsize(13)
        ax.yaxis.get_label().set_fontsize(13)
        ax.tick_params(axis='x', labelsize=12)
        ax.tick_params(axis='y', labelsize=12)
        plt.tight_layout()
        plt.show()

def show_synthdata_res(varying='data', eval_method='f1'):
    """
    Function that displays the PCC and SRCC results on the synthetic datasets.
    [WARNING] You must run the run_synthdata for all data variations or all ft variations and the same values for 
              classif and eval_method first!
    @varying: String, defines what synthetic datasets sjpimd be considered between "data" or features. 
              If == "data", results of the synthetic datasets with varying data distribution.
              If != "data", results of the synthetic datasets with varying feature components.
    @eval_method: String, defines the classification evaluation metric used : either "f1" for f1-score or gmean.
    """
    file_a = './table/results/corr_synth_varying_' + varying + '_a_' + str(eval_method) + '.csv'
    file_b = './table/results/corr_synth_varying_' + varying + '_b_' + str(eval_method) + '.csv'
    file_c = './table/results/corr_synth_varying_' + varying + '_c_' + str(eval_method) + '.csv'
    files = [file_a, file_b, file_c]
    
    for file in files:
        print(file)
        results = pd.read_csv(file)
        if varying == 'data':
            metrics = ["IR", "AdjIR", "C1", "C2", "ID", "LRID", "IF", "SIMBA"]
        else:
            metrics = ["AdjIR", "SIMBA"]
        for metric in metrics:
            all_pcc = results['PCC_' + metric].to_list()
            all_srcc = results['SRCC_' + metric].to_list()
            print('PCC')
            print('Correlation ' + eval_method + ' scores & ' + metric + ': ' + str(np.mean(all_pcc)) + ' +- ' + str(np.std(all_pcc)))
            print("SRCC")
            print('Correlation ' + eval_method + ' scores & ' + metric + ': '  + str(np.mean(all_srcc)) + ' +- ' + str(np.std(all_srcc)))
            input()

def data_complexity_analysis():
    """
    Function that runs the data complexity analysis for all imbalance measures.
    It includes 20 data complexity measures, calculated with the pymfe package.
    The results are saved in a table named 'data_complexity_analysis.csv' in the table/results/ folder.
    """
    l_datasets = ['abalone', 'cardio10', 'cardio3', 'contraceptive', 'dermatology', 'ecoli', 'glass', 
                  'hayes-roth', 'knowledge', 'landsat', 'lenses', 'lymphography', 'new-thyroid', 'obesity', 
                  'pageblocks', 'penbased', 'soybean', 'steel', 'student', 'theorem', 'thyroid', 'wallfollowing', 
                  'webphishing', 'wholesale', 'wine', 'wine-quality', 'yeast', 'abalone-20_vs_8-9-10', 
                  'abalone9-18', 'banknote', 'bankruptcy', 'breastcancer', 'cleveland-0_vs_4', 'credit', 
                  'dermatology-6', 'ecoli-0-1-4-7_vs_2-3-5-6', 'ecoli1', 'glass-0-1-5_vs_2', 'glass2', 'glass4', 
                  'glass5', 'glass6', 'ionosphere', 'led7digit-0-2-4-5-6-7-8-9_vs_1', 'new-thyroid1', 
                  'page-blocks0', 'poker-9_vs_7',  'segment0', 'spambase', 'spect-heart', 'vehicle0', 'vowel0', 
                  'winequality-red-3_vs_5', 'wisconsin', 'yeast-1-2-8-9_vs_7', 'yeast3', 'yeast4', 'yeast5',
                  'drybean', 'room', 'htru2', 'purchase', 'chess', 'connect-4', 'loc_build', 'loc_floor', 'adult',
                   'skin', 'shuttle']

    l_col_label = ["label"] * len(l_datasets)
    l_col_ID = [["id"]] * len(l_datasets)
    data_complexity_ft = ['cls_coef', 'density', 'f1', 'f1v', 'f2', 'f3', 'f4', 'hubs', 'l1', 'l2', 'l3', 'lsc', 
                          'n1', 'n2', 'n3', 'n4', 't1', 't2', 't3', 't4']
    data_complexity_ft_mean = ['f1', 'f1v', 'f2', 'f3', 'f4', 'hubs', 'l1', 'l2', 'l3', 'n2', 'n3', 'n4']
    df = pd.DataFrame(columns=['DS', 'IR', 'AdjIR', 'C1', 'C2', 'ID', 'LRID', 'IF', 'SIMBA', 'f1-score'] + data_complexity_ft)
    safe_file = pd.read_csv('./table/results/real_SVM_f1.csv')
    for i in range(len(l_datasets)):
        dataset = l_datasets[i]
        print("DATASET: " + str(dataset))
        col_label = l_col_label[i]
        col_ID = l_col_ID[i]
        ds_metrics_f1 = safe_file[safe_file['DS'] == dataset][['IR', 'AdjIR', 'C1', 'C2', 'ID', 'LRID', 'IF', 'SIMBA', 'f1']]
        path = "./data/" + dataset + "/"
        X = pd.read_csv(path + "df_allfeatures.csv") # features data
        y = pd.read_csv(path + "labels.csv") # labels data

        if dataset in ['connect-4', 'loc_build', 'loc_floor', 'adult', 'skin', 'shuttle']:
            sample_X = X.sample(frac=0.1, random_state=42).sort_index()
            sample_y = y.drop(index=[ind for ind in y.index if not (ind in sample_X.index)])
            X = sample_X
            y = sample_y            

        X = X.drop(columns=col_ID)
        y = y[col_label] 
        kf = StratifiedKFold(n_splits=5) # 5-fold CV
        tmp_df_5cv = pd.DataFrame(columns=data_complexity_ft)
        for k, (train_index, test_index) in enumerate(kf.split(X, y)):
            print('Fold ' + str(k+1))
            X_train = X.iloc[train_index]
            y_train = y.iloc[train_index]
            # with pymfe
            mfe = MFE(groups=['complexity'], features=data_complexity_ft)
            mfe.fit(np.asarray(X_train.astype(float)), np.asarray(y_train))
            ft = mfe.extract()
            ft_names = ft[0]
            ft_values = ft[1]
            row_tmp = {}
            for f in data_complexity_ft:
                if f in data_complexity_ft_mean:
                    name_ft = f + '.mean'
                else:
                    name_ft = f
                ft_val = ft_values[ft_names.index(name_ft)]
                row_tmp[f] = ft_val
            tmp_df_5cv.loc[len(tmp_df_5cv)] = row_tmp
        row = {'DS': dataset, 'IR': ds_metrics_f1['IR'].values[0], 'AdjIR': ds_metrics_f1['AdjIR'].values[0], 
               'C1': ds_metrics_f1['C1'].values[0], 'C2': ds_metrics_f1['C2'].values[0], 'ID': ds_metrics_f1['ID'].values[0], 
               'LRID': ds_metrics_f1['LRID'].values[0], 'IF': ds_metrics_f1['IF'].values[0], 
               'SIMBA': ds_metrics_f1['SIMBA'].values[0], 'f1-score': ds_metrics_f1['f1'].values[0]}
        for f in data_complexity_ft:
            all_val_5cv = tmp_df_5cv[f].to_list()
            all_val_5cv = [v for v in all_val_5cv if not np.isnan(v)]
            row[f] = np.mean(all_val_5cv)
        df.loc[len(df)] = row
    print(df)
    # input()
    df.to_csv('./table/results/data_complexity_analysis.csv', index=False)

def show_data_complexity():
    """
    Function that prints the results of the data complexity analysis.
    [WARNING] You must have run the data_complexity_analysis() function first for it to work.
    """
    file = './table/results/data_complexity_analysis.csv'
    metrics_imbalance = ['IR', 'AdjIR', 'C1', 'C2', 'ID', 'LRID', 'IF', 'SIMBA']
    metrics_complexity = ['f1', 'f1v', 'f2', 'f3', 'f4', # Feature-based measures
                          'l1', 'l2', 'l3',  # Linearity measures
                          'n1', 'n2', 'n3', 'n4', 't1', 'lsc', # Neighbourhood measures
                          'density', 'cls_coef', 'hubs', # Network measures
                          't2', 't3', 't4'] # Dimensionality measures
    eval_method = 'f1-score'

    for m_complex in metrics_complexity:
        print('TEST ' + m_complex)
        data_complexity = pd.read_csv(file)
        nan_val = data_complexity[m_complex].isna().sum()
        if nan_val != 0:
            data_complexity = data_complexity.drop(index=data_complexity[data_complexity[m_complex].isna()].index) # drop rows with NaN values for this data complexity measure
        print('NB DS: ' + str(len(data_complexity)))      
        data_m_complex = data_complexity[m_complex].to_list()
        data_eval = data_complexity[eval_method].to_list()
        data_m_complex = np.array(data_m_complex).reshape(-1, 1)
        data_eval  = np.array(data_eval).reshape(-1, 1)
        corr_perf = np.round(r_regression(data_m_complex, data_eval), 2)[0]
        p_val_perf = np.round(f_regression(data_m_complex, np.array(data_eval))[1][0],5)
        print('Correlation between '+ m_complex + ' and ' + eval_method + ': ' + str(corr_perf) + ', p-value: ' + str(p_val_perf))
        for m_imbalance in metrics_imbalance:
            data_m_imb = data_complexity[m_imbalance].to_list()
            data_m_imb = np.array(data_m_imb).reshape(-1, 1)
            corr_imb = np.round(r_regression(data_m_complex, data_m_imb), 2)[0]
            p_val_imb= np.round(f_regression(data_m_complex, np.array(data_m_imb))[1][0],5)
            print('Correlation between '+ m_complex + ' and ' + m_imbalance + ': ' + str(corr_imb) + ', p-value: ' + str(p_val_imb))
        input()


if __name__ == "__main__":
    warnings.simplefilter('ignore', category=ConstantInputWarning)
    warnings.simplefilter('ignore', category=RuntimeWarning)
    warnings.simplefilter('ignore', category=UserWarning)
    
    # Possible actions are:
    #   - table_ds: get descriptive tables of the real datasets that are used
    #   - synth_data: run study on synthetic datasets with varying data distribution
    #   - synth_ft: run study on synthetic datasets with varying feature components
    #   - synth_data_plot: plot results of action synth_data
    #   - synth_ft_plot: plot results of action synth_ft
    #   - real_data: run study on the 70 real datasets
    #   - real_data_plot: plot results of action real_data (regression lines for all, binary, and multiclass datasets)
    #   - ablation: run ablation study with SIMBA 4 versions: whole, no normalization, no feature importance, and no
    #                                                         feature redundancy
    #   - ablation_plot: plot results of action ablation
    #   - complexity: run the control study with data complexity measures

    action = 'table_ds'
    eval = "f1" # "f1" or "gmean"
    N = 100 # nb runs for synthetic datasets

    if action == "table_ds":
        only_class_freq = False # if False, generates table with only mean, sd for class frequency
                                # if True, generates table with detailed class frequencies
        # List of all datasets = separated between multiclass & binary datasets
        l_datasets = ['abalone', 'balance', 'cardio10', 'cardio3', 'chess', 
                    'connect-4', 'contraceptive', 'dermatology', 'drybean', 'ecoli', 
                    'glass', 'hayes-roth', 'knowledge', 'landsat', 'lenses', 
                    'loc_build', 'loc_floor', 'lymphography', 'new-thyroid', 'obesity', 
                    'pageblocks', 'penbased', 'room', 'shuttle', 'soybean', 
                    'steel', 'student', 'theorem', 'thyroid', 'wallfollowing', 
                    'webphishing', 'wholesale', 'wine', 'wine-quality', 'yeast']
        l_datasets_bin = ['abalone-20_vs_8-9-10', 'abalone9-18', 'adult', 'banknote', 'bankruptcy',
                        'breastcancer', 'cleveland-0_vs_4', 'credit', 'dermatology-6', 'ecoli-0-1-4-7_vs_2-3-5-6', 
                        'ecoli1', 'glass-0-1-5_vs_2', 'glass2', 'glass4', 'glass5', 
                        'glass6', 'htru2', 'ionosphere', 'led7digit-0-2-4-5-6-7-8-9_vs_1', 'new-thyroid1', 
                        'page-blocks0', 'poker-9_vs_7', 'purchase', 'segment0', 'skin',
                        'spambase', 'spect-heart', 'vehicle0', 'vowel0', 'winequality-red-3_vs_5', 
                        'wisconsin', 'yeast-1-2-8-9_vs_7', 'yeast3', 'yeast4', 'yeast5']
        all_ds = l_datasets + l_datasets_bin
        all_ds.sort()
        print(all_ds)
        make_DS_descriptive_table(all_ds, only_class_freq=only_class_freq)
    
    elif action == "synth_data":
        # SYNTHETIC DATASETS
        # VARYING DATA DISTRIBUTION, BUT FIXED FEATURES (10 FT, 2 INFORMATIVE)
        run_synthdata(N = N, data_distrib = [100, 5, 0], feature_distrib = [2, 0, 8], tot_ft = 10, i_range=[5, 101, 1], 
                      title='Datasets (a): (100, 5, i)', corr=True, varying="data", eval_method=eval)
        run_synthdata(N = N, data_distrib = [50, 50, 0], feature_distrib = [2, 0, 8], tot_ft = 10, i_range=[5, 101, 1], 
                      title='Datasets (b): (50, 50, i)', corr=True, varying="data", eval_method=eval)
        run_synthdata(N = N, data_distrib = [10, 50, 100], feature_distrib = [2, 0, 8], tot_ft = 10, i_range=[5, 101, 1], 
                      title='Datasets (c): (10i, 50i, 100i)', corr=True, varying="data", eval_method=eval)
        show_synthdata_res(varying='data', eval_method=eval)

    elif action == "synth_ft":
        # SYNTHETIC DATASETS
        # VARYING FEATURE DISTRIBUTION, BUT FIXED DATA DISTRIBUTION (400, 75, 25)
        run_synthdata(N = N, data_distrib = [400, 75, 25], feature_distrib = [-1, 0, 0], tot_ft = 50, i_range=[2, 51, 1], 
                      title='Datasets with i out of 50 informative features, dataset (400, 75, 25)', 
                      corr=True, varying="ft", eval_method=eval) 
        run_synthdata(N = N, data_distrib = [400, 75, 25], feature_distrib = [5, -1, 0], tot_ft = -1, i_range=[0, 46, 1], 
                      title='Datasets with 5+i features, 5 informative, and i redundant, dataset (400, 75, 25)', 
                      corr=True, varying="ft", eval_method=eval)
        run_synthdata(N = N, data_distrib = [400, 75, 25], feature_distrib = [0, -1, 0], tot_ft = 50, i_range=[0, 49, 1], 
                      title='Datasets with 50 features, 50-i informative, and i redundant, dataset (400, 75, 25)', 
                      corr=True, varying="ft", eval_method=eval)
        show_synthdata_res(varying='ft', eval_method=eval)
        
    elif action == "synth_data_plot":
        # PLOTS
        run_synthdata_plot(N=N, varying='data', eval_method=eval, from_file=True)
    
    elif action == "synth_ft_plot":
        # PLOTS
        run_synthdata_plot(N=N, varying='ft', eval_method=eval, from_file=True)

    elif action == "real_data":
        # REAL DATASETS with 5 common ML classifiers: SVM, LDA, RF, kNN, MLP
        clfs = ["SVM", "LDA", "RF", "kNN", "MLP"]
        for clf in clfs:
            print(clf + " running...")
            run_realdata(classif=clf, eval_method=eval)
            print(clf + " done!")
            print()

    elif action == "real_data_plot":
        # Regression plot real DS, MULTI+BIN, BIN, MULTI
        clfs = ["SVM", "LDA", "RF", "kNN", "MLP"]
        for clf in clfs:
            show_realdata_res(classif=clf, eval_method=eval)

    elif action == 'ablation':
        print("Abaltion study running...")
        run_real_ablation(classif = 'SVM', eval_method='f1')
        print("Done!")
    
    elif action == 'ablation_plot':
        plot_real_ablation(classif="SVM", eval_method="f1")
    
    elif action == "complexity":
        data_complexity_analysis()
        show_data_complexity()