features_analysis.py

# Copyright CNRS/Inria/UNS
# Contributor(s): Eric Debreuve (since 2019), Morgane Nadal (2020)
#
# eric.debreuve@cnrs.fr
#
# This software is governed by the CeCILL  license under French law and
# abiding by the rules of distribution of free software.  You can  use,
# modify and/ or redistribute the software under the terms of the CeCILL
# license as circulated by CEA, CNRS and INRIA at the following URL
# "http://www.cecill.info".
#
# As a counterpart to the access to the source code and  rights to copy,
# modify and redistribute granted by the license, users are provided only
# with a limited warranty  and the software's author,  the holder of the
# economic rights,  and the successive licensors  have only  limited
# liability.
#
# In this respect, the user's attention is drawn to the risks associated
# with loading,  using,  modifying and/or developing or reproducing the
# software by the user in light of its specific status of free software,
# that may mean  that it is complicated to manipulate,  and  that  also
# therefore means  that it is reserved for developers  and  experienced
# professionals having in-depth computer knowledge. Users are therefore
# encouraged to load and test the software's suitability as regards their
# requirements in conditions enabling the security of their systems and/or
# data to be ensured and,  more generally, to use and operate it in the
# same conditions as regards security.
#
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.

import pandas as pd_
import numpy as np_
import matplotlib.pyplot as pl_
from matplotlib.lines import Line2D
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn import tree
import scipy as si_
import seaborn as sb_
import itertools

import sys as sy_
from typing import List


def PCAOnDF(df: pd_.DataFrame(),
            target: str,
            targets: List[str],
            colors: List[str],
            save_name: str = None,
            title: str = "",
            plot_duration: bool = True,
            save_in: str = None,
            three_D: bool = False,
            ) -> list:
    '''
    Perform 2D or 3D PCA on the CHO-DIO dataframe.
    Save explained variance ratio as csv, plot and save the PCAs.
    '''
    # Separating the features from their conditions and durations
    all_target = pd_.DataFrame(df.loc[:, [target]].values, columns=[target])
    df_all = df.drop([target, "duration"], axis=1)

    # Standardize the data
    scaler = StandardScaler()
    scaler.fit(df_all)
    stand_df = scaler.transform(df_all)

    # Create the PCA and fit the data
    pca = PCA(n_components=2)
    principal_components = pca.fit_transform(stand_df)

    # Printing and saving the explained variance ratio
    print(f"PCA explained variance ratio ({save_name}): ", pca.explained_variance_ratio_)
    pca_exp_var_df = pd_.DataFrame(pca.explained_variance_ratio_, columns=["pca explained variance ration"])
    pca_exp_var_df = pca_exp_var_df.rename(index={0: "principal component 1", 1: "principal component 2"})
    pca_exp_var_df.to_csv(f"{save_in}\\pca_explained_variance_ratio_{save_name}.csv")

    # Creating the principal component df
    principal_df = pd_.DataFrame(data=principal_components, columns=['principal component 1', 'principal component 2'])
    # Give the final df containing the principal component and their condition
    final_df = pd_.concat([principal_df, all_target[target]], axis=1)

    # Plot 2D
    fig = pl_.figure(figsize=(8, 8))
    ax = fig.add_subplot(1, 1, 1)
    ax.set_xlabel(f'Principal Component 1 - ratio = {round(pca.explained_variance_ratio_[0],2)}', fontsize=15)
    ax.set_ylabel(f'Principal Component 2 - ratio = {round(pca.explained_variance_ratio_[1],2)}', fontsize=15)
    ax.set_title(f'2 component PCA{title}', fontsize=20)
    for tgt, color in zip(targets, colors):
        idx = final_df[target] == tgt
        ax.scatter(final_df.loc[idx, 'principal component 1']
                   , final_df.loc[idx, 'principal component 2']
                   , c=color
                   , s=30)
    ax.legend(targets)
    ax.grid()

    if save_name is not None:
        pl_.savefig(f"{save_in}\\PCA_{save_name}.png")

    if plot_duration:
        # Make separated plots for each duration of the experiments
        fig = pl_.figure(figsize=(8, 8))
        ax = fig.add_subplot(1, 1, 1)

        ax.set_xlabel(f'Principal Component 1 - ratio = {round(pca.explained_variance_ratio_[0],2)}', fontsize=15)
        ax.set_ylabel(f'Principal Component 2 - ratio = {round(pca.explained_variance_ratio_[1],2)}', fontsize=15)
        ax.set_title(f'2 component PCA{title}', fontsize=20)

        new_tgts = [["CHO", "1H"], ["CHO", "3H"], ["CHO", "6H"], ["DIO", "1H"], ["DIO", "3H"], ["DIO", "6H"]]
        final_df = pd_.concat([df[["condition", "duration"]], principal_df], axis=1)

        for tgt, color in zip(new_tgts, ["lavender", "royalblue", "navy", "lightcoral", "firebrick", "red"]):
            new_df = final_df.loc[(final_df["condition"] == tgt[0]) & (final_df["duration"] == tgt[1])]
            new_df = new_df.drop(["condition", "duration"], axis=1)

            ax.scatter(new_df['principal component 1']
                       , new_df['principal component 2']
                       , c=color
                       , s=30)

        ax.legend(new_tgts)
        ax.grid()

        if save_name is not None:
            pl_.savefig(f"{save_in}\\PCA_duration_{save_name}.png")

    if three_D:
        # Create the 3D PCA and fit the data
        pca3d = PCA(n_components=3)
        principal_components3d = pca3d.fit_transform(stand_df)

        # Print and Save the explained variance ratio
        print(f"3D PCA explained variance ratio ({save_name}): ", pca3d.explained_variance_ratio_)
        pca3d_exp_var_df = pd_.DataFrame(pca3d.explained_variance_ratio_, columns=["3d pca explained variance ration"])
        pca3d_exp_var_df = pca3d_exp_var_df.rename(index={0: "principal component 1", 1: "principal component 2", 2: "principal component 3"})
        pca3d_exp_var_df.to_csv(f"{save_in}\\three_d_pca_explained_variance_ratio_{save_name}.csv")

        # Creating the principal component df
        principal3d_df = pd_.DataFrame(data=principal_components3d,
                                       columns=['principal component 1', 'principal component 2', 'principal component 3'])
        # Give the final df containing the principal component and their condition
        final3d_df = pd_.concat([principal3d_df, all_target[target]], axis=1)

        # Plot
        fig = pl_.figure(figsize=(8, 8))
        ax = fig.add_subplot(111, projection='3d')

        ax.set_xlabel(f'Principal Component 1 - ratio = {round(pca3d.explained_variance_ratio_[0],2)}', fontsize=15)
        ax.set_ylabel(f'Principal Component 2 - ratio = {round(pca3d.explained_variance_ratio_[1],2)}', fontsize=15)
        ax.set_zlabel(f'Principal Component 3 - ratio = {round(pca3d.explained_variance_ratio_[2],2)}', fontsize=15)
        ax.set_title(f'3 component PCA{title}', fontsize=20)

        for tgt, color in zip(targets, colors):
            idx = final3d_df[target] == tgt

            ax.scatter(final3d_df.loc[idx, 'principal component 1']
                       , final3d_df.loc[idx, 'principal component 2']
                       , final3d_df.loc[idx, 'principal component 3']
                       , c=color
                       , s=30)

        ax.legend(targets)
        ax.grid()

        if save_name is not None:
            pl_.savefig(f"{save_in}\\three_d_PCA_{save_name}.png")

        if plot_duration:
            fig = pl_.figure(figsize=(8, 8))
            ax = fig.add_subplot(111, projection="3d")

            ax.set_xlabel(f'Principal Component 1 - ratio = {round(pca3d.explained_variance_ratio_[0], 2)}', fontsize=15)
            ax.set_ylabel(f'Principal Component 2 - ratio = {round(pca3d.explained_variance_ratio_[1],2)}', fontsize=15)
            ax.set_zlabel(f'Principal Component 3 - ratio = {round(pca3d.explained_variance_ratio_[2],2)}', fontsize=15)
            ax.set_title(f'3 component PCA{title}', fontsize=20)

            final3d_df = pd_.concat([df[["condition", "duration"]], principal3d_df], axis=1)

            for tgt, color in zip(new_tgts, ["lavender", "royalblue", "navy", "lightcoral", "firebrick", "red"]):
                new3d_df = final3d_df.loc[(final_df["condition"] == tgt[0]) & (final3d_df["duration"] == tgt[1])]
                new3d_df = new3d_df.drop(["condition", "duration"], axis=1)

                ax.scatter(new3d_df['principal component 1']
                           , new3d_df['principal component 2']
                           , new3d_df['principal component 3']
                           , c=color
                           , s=30)

            ax.legend(new_tgts)
            ax.grid()

            if save_name is not None:
                pl_.savefig(f"{save_in}\\three_d_PCA_duration_{save_name}.png")


def KmeansOnDF(df: pd_.DataFrame(),
               nb_clusters: tuple,
               target: str,
               plot_bar: bool = True,
               rep_on_image: bool = False,
               labeled_somas=None,
               elbow: bool = False,
               intracluster_var: bool = True,
               save_name: str = None,
               save_in: str = None,
               title: str = "",
               duration: bool = False,
               features_distribution: bool = True,
               ) -> KMeans:
    '''
    Perform kmeans on the pandas dataframe. Can find the best number of cluster with elbow method,
    find the intracluster variance, and represent the result on the initial images.
    Plot barplots with percentage of each label for each condition.
    '''
    # Separating the features from their conditions and durations
    all_target = pd_.DataFrame(df.loc[:, [target]].values, columns=[target])
    df2 = df.drop([target, "duration"], axis=1)
    df_scalar = df.drop(["duration", "condition"], axis=1)

    # Data standardization
    scaler = StandardScaler()
    stand_df = scaler.fit_transform(df2)

    # Best number of clusters using Elbow method
    if elbow:
        wcss = []  # within cluster sum of errors(wcss)
        for i in range(1, 24):
            kmeans = KMeans(n_clusters=i, init='k-means++', max_iter=300, n_init=10, random_state=0)
            kmeans.fit(stand_df)
            wcss.append(kmeans.inertia_)
        pl_.figure()
        pl_.plot(range(1, 24), wcss)
        pl_.plot(range(1, 24), wcss, 'bo')
        pl_.title('Elbow Method')
        pl_.xlabel('Number of clusters')
        pl_.ylabel('WCSS')
        pl_.show(block=True)
        pl_.close()

    # Kmeans with x clusters
    for nb_cluster in nb_clusters:
        kmeans = KMeans(n_clusters=nb_cluster, init='k-means++', max_iter=300, n_init=10, random_state=0)
        kmeans.fit_predict(stand_df)
        label_df = pd_.DataFrame(data=kmeans.labels_, columns=['label'])
        lab_cond_df = pd_.concat([label_df, all_target[target]], axis=1)

        # Intracluster variance
        if intracluster_var:
            var_df = pd_.DataFrame(stand_df)
            var = IntraClusterVariance(var_df, kmeans, nb_cluster, save_name)
            var_df = pd_.DataFrame(var, columns=["intracluster variance"])
            var_df = var_df.rename({cluster: f"label {cluster}" for cluster in range(nb_cluster)})
            if save_in:
                var_df.to_csv(f"{save_in}\\intracluster_variance_kmeans_{save_name}_k={nb_cluster}.csv")

        # TODO Intersection of ellipsoids of the covariance matrices
        #  - ellipsoid equation: (x-mean)(1/Mcov)(x-mean).T <= 1

        # Barplot
        if plot_bar:
            if duration:
                fig = pl_.figure(figsize=(16, 8))
                ax = fig.add_subplot(1, 1, 1)
                (lab_cond_df.groupby("condition")["label"].value_counts(normalize=True).mul(100).rename(
                    "Percentage (%)").reset_index().pipe((sb_.barplot, "data"), x="condition", y="Percentage (%)",
                                                         hue="label", palette=sb_.color_palette("deep", n_colors=nb_cluster)))
                ax.set_ylim(0, 100)
                ax.set_title(f'Distribution of the clustering labels according to conditions{title}_k={nb_cluster}', fontsize=11)
                ax.grid()
                if save_name is not None:
                    pl_.savefig(f"{save_in}\\Hist_Clustering_{save_name}_k={nb_cluster}.png")
                    # pl_.show(block=True)
                    # pl_.close()
            else:
                ## Batthacharyya similarity
                groupbycond = lab_cond_df.groupby("condition")
                hist = []
                for cond, values in groupbycond:
                    hist.append(np_.histogram(values['label'], bins=2))
                bhatt = BhattacharyyaSimilarity(hist[0][0], hist[1][0])
                print("Bhattacharyya similarity btw CHO and DIO, all times: ", bhatt)

                ## Plot the barplots
                fig = pl_.figure(figsize=(16, 8))
                ax= fig.add_subplot(1, 1, 1)
                (groupbycond["label"].value_counts(normalize=True).mul(100).rename(
                    "Percentage (%)").reset_index().pipe((sb_.barplot, "data"), x="condition", y="Percentage (%)",
                                                         hue="label", palette=sb_.color_palette("deep", n_colors=nb_cluster)))
                ax.set_ylim(0, 100)
                ax.grid()
                ax.set_title(f'Distribution of the clustering labels according to conditions{title}_k={nb_cluster}', fontsize=11)
                if save_name is not None:
                    pl_.savefig(f"{save_in}\\Hist_Clustering_{save_name}_k={nb_cluster}.png")
                    # pl_.show(block=True)
                    # pl_.close()
                lab_cond_df2 = pd_.concat((lab_cond_df, df[["duration"]]), axis=1)

                bhattacharyya = []
                for dur in df["duration"].unique():
                    dur_df = lab_cond_df2.loc[lab_cond_df2["duration"] == dur]

                    ## Bhattacharyya
                    groupbycond = dur_df.groupby("condition")
                    hist = []
                    for cond, values in groupbycond:
                        hist.append(np_.histogram(values['label'], bins=2))
                    bhatt = BhattacharyyaSimilarity(hist[0][0], hist[1][0])
                    print(f"Bhattacharyya similarity btw CHO and DIO {dur}: ", bhatt)
                    bhattacharyya.append(bhatt)

                    ## Barplot
                    fig = pl_.figure(figsize=(16, 8))
                    ax = fig.add_subplot(1, 1, 1)
                    to_plot = dur_df.groupby("condition")["label"].value_counts(normalize=True).mul(100).rename("Percentage (%)").reset_index()
                    n_labels = to_plot["label"].unique()
                    n_labels.sort()
                    colors = sb_.color_palette("deep", n_colors=nb_cluster)
                    color_to_use = list(colors[i] for i in n_labels)
                    (to_plot.pipe((sb_.barplot, "data"), x="condition", y="Percentage (%)", hue="label", palette=color_to_use))
                    ax.set_ylim(0, 100)
                    ax.grid()
                    ax.set_title(f'Distribution of the clustering labels according to conditions{title}_k={nb_cluster} - duration = {dur}', fontsize=11)
                    if save_name is not None:
                        pl_.savefig(f"{save_in}\\Hist_Clustering_{save_name}_dur_{dur}_k={nb_cluster}.png")
                        # pl_.show(block=True)
                        pl_.close()

                if bhattacharyya is not None:
                    fig = pl_.figure(figsize=(16, 8))
                    ax = fig.add_subplot(1, 1, 1)
                    ax.plot([f"{duration}" for duration in df["duration"].unique()],
                            [bhattacharyya[i] for i in range(3)], "ko")
                    ax.plot([f"{duration}" for duration in df["duration"].unique()],
                            [bhattacharyya[i] for i in range(3)], "g")
                    ax.set_xlabel("Duration")
                    ax.set_yscale("log")
                    ax.set_ylabel("Bhattacharyya distance")
                    ax.grid()
                    ax.set_title(f'Bhattacharyya distance between CHO and DIO{title}_k={nb_cluster}',
                                 fontsize=11)

                    if save_name is not None:
                        pl_.savefig(f"{save_in}\\Bhattacharyya_distance_{save_name}_k={nb_cluster}.png")
                        pl_.close()

                if save_name is not None:
                    pl_.savefig(f"{save_in}\\Hist_duration_Clustering_{save_name}_k={nb_cluster}.png")
                    print(f"Saved in {save_in}")
                    # pl_.show(block=True)
                    pl_.close()

        if features_distribution:
            for column in df_scalar.columns:
                cond_col_df = pd_.concat((df_scalar[column], lab_cond_df), axis=1)

                conditions = ['CHO', 'DIO']
                labels = np_.arange(nb_cluster)

                fig = pl_.figure(constrained_layout=False)
                gs = fig.add_gridspec(ncols=1, nrows=1)
                ax1 = fig.add_subplot(gs[0, 0])

                # Plot a histogram and kernel density estimate
                # print(f"Plot a histogram and kernel density estimate for feature {column}")
                for comb, color in zip(itertools.product(conditions, labels), sb_.color_palette("deep", n_colors=len(labels)*len(conditions))):
                    to_plot = cond_col_df.loc[(cond_col_df['condition'] == comb[0])]
                    to_plot = to_plot.loc[(to_plot['label'] == comb[1])]

                    # Kernel estimate of the histogram
                    sb_.distplot(to_plot[[column]], hist= False, color=color, ax=ax1)

                lines = [Line2D([0], [0], color=c, linewidth=3, linestyle='-') for c in sb_.color_palette("deep", n_colors=len(labels)*len(conditions))]
                lb = list(itertools.product(conditions, labels))
                ax1.legend(lines, lb)
                ax1.set_title(f'{column} distribution{title}', fontsize=11)

                pl_.tight_layout()

                if save_in is not None:
                    pl_.savefig(f"{save_in}\\feat_kernel_estimate_{column}_k={nb_cluster}.png")
                    pl_.close()

                # Boxplot
                for comb, color in zip(itertools.product(conditions, labels), sb_.color_palette("deep", n_colors=len(labels)*len(conditions))):
                    fig = pl_.figure(constrained_layout=False, figsize=(10, 3))
                    gs = fig.add_gridspec(ncols=1, nrows=1)
                    ax1 = fig.add_subplot(gs[0, 0])

                    to_plot = cond_col_df.loc[(cond_col_df['condition'] == comb[0])]
                    to_plot = to_plot.loc[(to_plot['label'] == comb[1])]

                    sb_.boxplot(to_plot[[column]], color=color, ax=ax1)

                    ax1.set_xlim(min(cond_col_df[column]), max(cond_col_df[column]))
                    lines = [Line2D([0], [0], color=color, linewidth=3, linestyle='-')]
                    ax1.legend(lines, [comb])
                    ax1.set_title(f'{column} boxplot{title}', fontsize=11)

                    pl_.tight_layout()

                    if save_in is not None:
                        pl_.savefig(f"{save_in}\\feat_boxplot_{column}_{comb}_k={nb_cluster}.png")
                        pl_.close()

                # Do the same thing but separating durations
                cond_col_df_dur = pd_.concat((cond_col_df, df["duration"]), axis=1)
                groupby_dur_ = cond_col_df_dur.groupby("duration")

                fig = pl_.figure(constrained_layout=False, figsize=(15, 10))
                gs = fig.add_gridspec(ncols=3, nrows=1)
                ax1 = fig.add_subplot(gs[0, 0])
                ax2 = fig.add_subplot(gs[0, 1])
                ax3 = fig.add_subplot(gs[0, 2])
                ax = [ax1, ax2, ax3]

                # Plot a histogram and kernel density estimate
                # print(f"Plot a histogram and kernel density estimate for feature {column}")
                x = 0
                for dur_, val in groupby_dur_:
                    for comb, color in zip(itertools.product(conditions, labels), sb_.color_palette("deep", n_colors=len(labels)*len(conditions))):
                        to_plot = val.loc[(val['condition'] == comb[0])]
                        to_plot = to_plot.loc[(to_plot['label'] == comb[1])]

                        # Kernel estimate of the histogram
                        sb_.distplot(to_plot[[column]], hist=False, color=color, ax=ax[x])

                    lines = [Line2D([0], [0], color=c, linewidth=3, linestyle='-') for c in sb_.color_palette("deep", n_colors=len(labels)*len(conditions))]
                    lb = list(itertools.product(conditions, labels))
                    ax[x].legend(lines, lb)
                    ax[x].set_title(f'{dur_}', fontsize=11)
                    # ax[x].set_xlim(min(val[column]), max(val[column]))
                    # ax[x].set_ylim()
                    ax[x].set_xlabel("Distribution kernel estimate")
                    ax[x].set_ylabel("Features values")
                    x += 1

                fig.suptitle(f'{column} distribution{title}')

                # pl_.tight_layout()

                if save_in is not None:
                    pl_.savefig(f"{save_in}\\feat_kernel_estimate_{column}_k={nb_cluster}_dur.png")
                    pl_.close()


        # Representation on the image
        if rep_on_image:
            RepresentationOnImages(labeled_somas, kmeans, nb_cluster)

    return kmeans


def IntraClusterVariance(df: pd_.DataFrame(), kmeans: KMeans(), nb_cluster: int, save_name: str = "") -> list:
    '''
    Return the intracluster variance of a given cluster found by kmeans.
    '''
    var = []
    # TODO change to .inertia_
    for cluster in range(nb_cluster):
        soma_cluster = [indx for indx, value in enumerate(kmeans.labels_) if value == cluster]
        mean_cluster = np_.average([df.iloc[row, :] for row in soma_cluster], axis=0)  # TODO .cluster_centers_
        variance = sum([np_.linalg.norm(df.iloc[row, :] - mean_cluster) ** 2 for row in soma_cluster]) / (len(soma_cluster) - 1)
        var.append(variance)

    print(f"Intracluster variance for {nb_cluster} clusters ({save_name}) :", var)

    return var


def RepresentationOnImages(labeled_somas, kmeans, nb_cluster):
    '''
    Represent the result of kmeans on labeled image. IN DVPT. Only available for a kmean intra-image.
    '''
    clustered_somas = labeled_somas.copy()
    clustered_somas = np_.amax(clustered_somas, axis=0)
    for indx, value in enumerate(kmeans.labels_):
        for indx_axe, axe in enumerate(clustered_somas):
            for indx_pixel, pixel in enumerate(axe):
                if pixel == indx + 1:
                    clustered_somas[indx_axe][indx_pixel] = value + 1
    pl_.imshow(clustered_somas, cmap="tab20")
    pl_.title(f"n cluster = {nb_cluster}")
    pl_.show(block=True)
    pl_.close()


def FeaturesStatistics(df: pd_.DataFrame(),
                       save_in: str = None,
                       title: str = "",
                       describe: bool = True,
                       heatmap: bool = True,
                       drop_feat: bool = True,
                       distribution: bool = True,
                       stat_test: bool = True,
                       decision_tree: bool = True,
                       ):
    '''
    Return the statistics allowing the user to choose the most relevant features to feed ML algorithms for ex.
    Statistical description, correlation heatmap, dropping correlated features or features with little difference between the two conditions,
    Plotting features distribution and boxplots, performing Kolmogorov-Smirnov two-sample test and Wilcoxon signed-rank two-sample test.
    Can draw a decision tree.
    '''
    #
    # Overview of the basic stats on each columns of df
    if describe:
        description = df.describe()
        if save_in is not None:
            description.to_csv(f"{save_in}\\df_stat_description.csv")

    # Data formatting
    df_scalar = df.drop(["duration", "condition"], axis=1)
    df_cond = df.drop(["duration"], axis=1)
    df_groupby_dur = df.groupby("duration")

    # Compute the correlation matrix
    corr_matrix = df_scalar.corr().abs()

    if heatmap:
        # Plot heat map with correlation matrix btw features
        print("Heat map with correlation matrix btw features")
        # Generate a mask for the upper triangle
        mask = np_.triu(np_.ones_like(corr_matrix, dtype=np_.bool))
        # Set up the figure
        fig, ax = pl_.subplots(figsize=(13, 9))
        # Generate a custom diverging colormap
        cmap = sb_.diverging_palette(220, 10, as_cmap=True)
        # Draw the heatmap with the mask and correct aspect ratio
        sb_.heatmap(corr_matrix, mask=mask, cmap=cmap, center=0, square=True, linewidths=.5, cbar_kws={"shrink": .5}, xticklabels=False)
        ax.set_title(f'Features correlation heat map{title}', fontsize=20)
        if save_in is not None:
            pl_.savefig(f"{save_in}\\Features correlation heat map{title}.png")
            pl_.close()

    if stat_test:
        # Create dictionaries to store the statistics and p-values
        dict_ks = {}
        dict_wx = {}
        dict_ks_dur = {}
        dict_wx_dur = {}

        for column in df_scalar.columns:
            # For each feature, perform a statistical test to know if the feature distribution or median is different btw CHO and DIO conditions
            cond_col_df = pd_.concat((df_scalar[column], df_cond["condition"]), axis=1)

            # Separate the conditions
            CHO = cond_col_df.loc[cond_col_df['condition'] == "CHO"]
            CHO = np_.asarray(CHO[column])
            DIO = cond_col_df.loc[cond_col_df['condition'] == "DIO"]
            DIO = np_.asarray(DIO[column])

            # Compare distribution between conditions (goodness of fit)
            ks = si_.stats.ks_2samp(CHO, DIO)
            dict_ks[column] = ks

            # Compare median between conditions
            wx = si_.stats.ranksums(CHO, DIO)
            dict_wx[column] = wx

            # Plot the p-values
            fig = pl_.figure(constrained_layout=False, figsize=(15,8))
            gs = fig.add_gridspec(ncols=2, nrows=1)
            ax1 = fig.add_subplot(gs[0, 0])
            ax2 = fig.add_subplot(gs[0, 1])
            #
            ax1.set_title(f"Kolmogorov-Smirnov p-value btw conditions - {column}", fontsize=13, pad=20)
            ax1.grid()
            ax1.set_xlabel("Duration")
            ax1.set_yscale("log")
            ax1.set_ylabel("p-value")
            #
            ax2.set_title(f"Wilcoxon signed-rank p-value btw conditions - {column}", fontsize=13, pad=20)
            ax2.grid()
            ax2.set_xlabel("Duration")
            ax2.set_yscale("log")
            ax2.set_ylabel("p-value")

            # Do the same thing than above but separate durations between conditions
            for duration, values in df_groupby_dur:
                if duration not in dict_ks_dur:
                    dict_ks_dur[duration] = {}
                    dict_wx_dur[duration] = {}

                duration_df = values.drop(["duration"], axis=1)
                CHO_ = duration_df.loc[duration_df['condition'] == "CHO"]
                DIO_ = duration_df.loc[duration_df['condition'] == "DIO"]
                CHO_ = CHO_.drop(["condition"], axis=1)
                DIO_ = DIO_.drop(["condition"], axis=1)
                CHO_ = np_.asarray(CHO_[column])
                DIO_ = np_.asarray(DIO_[column])

                # Compare distribution
                ks2 = si_.stats.ks_2samp(CHO_, DIO_)
                dict_ks_dur[duration][column] = ks2

                # Compare median
                wx2 = si_.stats.ranksums(CHO_, DIO_)
                dict_wx_dur[duration][column] = wx2

            ax1.plot([f"{duration}" for duration, _ in df_groupby_dur], [dict_ks_dur[duration][column][1] for duration, _ in df_groupby_dur], "g")
            ax1.plot([f"{duration}" for duration, _ in df_groupby_dur], [dict_ks_dur[duration][column][1] for duration, _ in df_groupby_dur], 'ko')
            ax2.plot([f"{duration}" for duration, _ in df_groupby_dur], [dict_wx_dur[duration][column][1] for duration, _ in df_groupby_dur], "y")
            ax2.plot([f"{duration}" for duration, _ in df_groupby_dur], [dict_wx_dur[duration][column][1] for duration, _ in df_groupby_dur], 'ko')

            pl_.tight_layout()

            if save_in:
                pl_.savefig(f"{save_in}\\p_values_duration_{column}")
                pl_.close()

        # Reformat the data to save it as csv
        df_ks = pd_.DataFrame.from_dict(data=dict_ks)
        df_ks = df_ks.rename(index={0: "Kolmogorov-Smirnov statistic", 1: "Kolmogorov-Smirnov p-value"})
        df_wx = pd_.DataFrame.from_dict(data=dict_wx)
        df_wx = df_wx.rename(index={0: "Wilcoxon statistic", 1: "Wilcoxon p-value"})

        df_ks_dur = pd_.DataFrame()
        df_wx_dur = pd_.DataFrame()
        for key in dict_ks_dur.keys():
            df_ks_dur = pd_.concat((df_ks_dur, pd_.DataFrame.from_dict(data=dict_ks_dur[key])))
            df_ks_dur = df_ks_dur.rename(index={0: f"{key} - Kolmogorov-Smirnov statistic", 1: f"{key} - Kolmogorov-Smirnov p-value"})
            df_wx_dur = pd_.concat((df_wx_dur, pd_.DataFrame.from_dict(data=dict_wx_dur[key])))
            df_wx_dur = df_wx_dur.rename(index={0: f"{key} - Wilcoxon statistic", 1: f"{key} - Wilcoxon p-value"})

        stat_tests_df1 = pd_.concat((df_ks, df_wx))  # KS and Wx all durations
        stat_tests_df2 = pd_.concat((df_ks_dur, df_wx_dur))  # KS and Wx for each duration

        if save_in:
            stat_tests_df1.to_csv(f"{save_in}\\stat_test_KS_wilcoxon_all.csv")
            stat_tests_df2.to_csv(f"{save_in}\\stat_test_KS_wilcoxon_for_each_duration.csv")

    if drop_feat:
        # Drop highly correlated features
        print("Drop highly correlated features")
        # Select upper triangle of correlation matrix
        upper = corr_matrix.where(np_.triu(np_.ones(corr_matrix.shape), k=1).astype(np_.bool))
        # Find index of feature columns with correlation greater than 0.9
        to_drop = [column for column in upper.columns if any(upper[column] > 0.9)]
        # Drop features
        drop_HCF_df = df_scalar.drop(df[to_drop], axis=1)

        # Drop column with null variance
        drop_HCF_df = drop_HCF_df.loc[:, drop_HCF_df.var() != 0.0]

        if stat_test:
            # drop non significant features (distribution and mean not different btw CHO and DIO)
            drop_HCF_df_6H = drop_HCF_df.copy()
            # based on all durations;
            to_drop =[column for column in drop_HCF_df.columns if (stat_tests_df1.loc["Kolmogorov-Smirnov p-value", column] > 1.0e-2) and (stat_tests_df1.loc["Wilcoxon p-value", column] > 1.0e-2)]
            drop_HCF_df = drop_HCF_df.drop(drop_HCF_df[to_drop], axis=1)

            # only based on 6H duration;
            to_drop_6H =[column for column in drop_HCF_df_6H.columns if (stat_tests_df2.loc["6H - Kolmogorov-Smirnov p-value", column] > 1.0e-3) and (stat_tests_df2.loc["6H - Wilcoxon p-value", column] > 1.0e-3)]
            drop_HCF_df_6H = drop_HCF_df_6H.drop(drop_HCF_df_6H[to_drop_6H], axis=1)

            drop_HCF_df_6H = pd_.concat((df[["condition", "duration"]], drop_HCF_df_6H), axis=1)

            if save_in:
                drop_HCF_df_6H.to_csv(f"{save_in}\\df_drop_highly_corr_feat_6H.csv")
                print(f"Selection of features with distribution and median of CHO vs DIO different in: {save_in}")

        # Add the condition and duration
        drop_HCF_df = pd_.concat((df[["condition", "duration"]], drop_HCF_df), axis=1)

        if save_in:
            drop_HCF_df.to_csv(f"{save_in}\\df_drop_highly_corr_feat.csv")
            print(f"Selection of low correlated features in: {save_in}")

    # Statistics for each features
    if distribution:
        for column in df_scalar.columns:
            cond_col_df = pd_.concat((df_scalar[column], df_cond["condition"]), axis=1)

            fig = pl_.figure(constrained_layout=False)
            gs = fig.add_gridspec(ncols=2, nrows=1)
            ax1 = fig.add_subplot(gs[0, 0])
            ax2 = fig.add_subplot(gs[0, 1])

            # Plot a histogram and kernel density estimate
            print(f"Plot a histogram and kernel density estimate for feature {column}")
            CHO = cond_col_df.loc[cond_col_df['condition'] == "CHO"]
            DIO = cond_col_df.loc[cond_col_df['condition'] == "DIO"]

            sb_.distplot(CHO[[column]], color="b", ax=ax1)
            sb_.distplot(DIO[[column]], color="r", ax=ax1)

            # Draw a boxplot
            print(f"Plot a boxplot for feature {column}")
            sb_.boxplot(data=df_cond, x="condition", y=column, hue="condition", palette=["b", "r"], ax=ax2)
            ax1.set_title(f'{column} distribution{title}', fontsize=11)
            ax2.set_title(f'{column} boxplot{title}', fontsize=11)

            pl_.tight_layout()

            if save_in is not None:
                pl_.savefig(f"{save_in}\\feat_distrib_{column}.png")
                pl_.close()

    # Decision tree
    if decision_tree:
        # Test btw CHO and DIO
        dt_df = df.drop(["condition", "duration"], axis=1)
        clf = tree.DecisionTreeClassifier(max_depth=4)
        clf = clf.fit(dt_df, df["condition"])
        fig = pl_.figure(figsize=(150, 65))
        tree.plot_tree(clf, feature_names=df.columns, class_names=["CHO", "DIO"], filled=True, rounded=True, fontsize=60)
        fig.suptitle(f"Decision tree all durations", fontsize=120)
        if save_in:
            pl_.savefig(f"{save_in}\\Decision_tree_all_durations_{title}.png")
            pl_.close()

        # Test btw CHO and DIO depending on duration
        for duration, values in df_groupby_dur:
            duration_df = values.drop(["duration", "condition"], axis=1)
            clf = tree.DecisionTreeClassifier(max_depth=3)
            clf = clf.fit(duration_df, values["condition"])

            fig = pl_.figure(figsize=(30, 16))
            tree.plot_tree(clf, feature_names=df.columns, class_names=["CHO", "DIO"], filled=True, rounded=True, fontsize=8)
            fig.suptitle(f"Decision tree {duration}", fontsize=16)
            if save_in:
                pl_.savefig(f"{save_in}\\Decision_tree_{duration}_{title}.png")
                pl_.close()


def BhattacharyyaSimilarity(h1, h2):
    return - np_.log(np_.sum(np_.sqrt(np_.multiply(Normalize(h1), Normalize(h2)))))


def Normalize(h):
    return h / np_.sum(h)


if __name__ == "__main__":
    # TODO: clean, reduce and optimize the code (many duplicates)
    #
    # os.chdir("path")

    ## If need to concatenate files:
    # all_filenames = [i for i in glob.glob('*.{}'.format("csv"))]
    # print(all_filenames)
    # df = pd_.concat([pd_.read_csv(f, index_col=0) for f in all_filenames])
    # df.to_csv(".\combined_features.csv")

    ## If use labeled somas:
    # labeled_somas = np_.load("path.npy")
    # df = pd_.read_csv(".\combined_features.csv", index_col=0)

    ## Parameters
    path = sy_.argv[1]
    save_in = sy_.argv[2]

    ## DF cleaning
    df0 = pd_.read_csv(f"{path}\\features.csv",
                      # index_col=0,
                      )
    df = df0.drop(["Unnamed: 0"], axis=1)

    ## Statistical analysis
    # For the moment drop the columns with non scalar values, and un-useful values
    # - TO BE CHANGED TODO (use distance metrics such as bhattacharyya coef, etc)
    df = df.drop(["soma uid",
                  "spherical_angles_eva", "spherical_angles_evb",
                  "hist_lengths", "hist_lengths_P", "hist_lengths_S",
                  "hist_curvature", "hist_curvature_P", "hist_curvature_S"],
                 axis=1)
    df = df.dropna(axis=0, how="any")

    # -- PCA with all the features
    print("\nALL FEATURES\n")
    # Between the two conditions, regardless the duration of experiment (2 conditions, all durations)
    PCAOnDF(df,
            target="condition",
            targets=["CHO", "DIO"],
            colors=["b", "r"],
            save_name="all_features",
            save_in=save_in,
            three_D=True,
            )

    # # Between the two conditions, for each duration (2 conditions, 3 durations)
    # groupby_duration = df.groupby("duration")
    # for duration, values in groupby_duration:
    #     ## duration: str, values: pd_.DataFrame()
    #     PCAOnDF(values,
    #             target="condition",
    #             targets=["CHO", "DIO"],
    #             colors=["b", "r"],
    #             save_name=f"{duration}_features",
    #             save_in=save_in,
    #             title=f" - {duration} Sample",
    #             plot_duration=False,
    #             three_D=True,
    #             )

    # -- K-means with all the features (2 conditions)

    # Test for multiple glial populations
    # Between the two conditions, regardless the duration of experiment (2 conditions, all durations)
    kmeans = KmeansOnDF(df,
                        target="condition",
                        nb_clusters=(2, 3, 4, 5),
                        elbow=True,
                        intracluster_var=True,
                        plot_bar=True,
                        save_name="all_features_multiple_pop",
                        save_in=save_in,
                        features_distribution=False,
                        )

    # # Between the two conditions, for each duration (2 conditions, 3 durations)
    # groupby_duration = df.groupby("duration")
    # for duration, values in groupby_duration:
    #     kmeans = KmeansOnDF(values,
    #                         target="condition",
    #                         nb_clusters=(2, 3, 4, 5),
    #                         elbow=False,
    #                         intracluster_var=True,
    #                         plot_bar=True,
    #                         save_name=f"{duration}_features_multiple_pop",
    #                         title=f" - {duration} Sample",
    #                         duration=True,
    #                         save_in=save_in,
    #                         )

    # -- Various plots to analyse the data and find discriminant features by statistical analysis
    print("\nFEATURE SELECTION\n")
    FeaturesStatistics(df,
                       save_in=save_in,
                       describe=True,
                       heatmap=True,
                       distribution=True,
                       stat_test=True,
                       drop_feat=True,
                       decision_tree=False,
                       )
    ## TODO: Enter selected features here
    # selected_features = []
    # selected_df = df[selected_features]
    ## TODO Or use the csv with dropped features
    try:
        selected_df = pd_.read_csv(f"{save_in}\\df_drop_highly_corr_feat.csv")
        # selected_df = pd_.read_csv(f"{save_in}\\df_drop_highly_corr_feat_6H.csv")
    except:
        raise("Only run the part until FeaturesStatistics included to generate df_drop_highly_corr_feat.csv, and then run the last part.")
        ## If an error raises, only run the part until FeaturesStatistics included, and then run the last part.

    # if other columns need to be dropped:
    try:
        to_drop = ["Unnamed: 0", "min_curvature"]
        selected_df = selected_df.drop(to_drop, axis=1)
    except:
        selected_df = selected_df.drop(["Unnamed: 0"], axis=1)

    # -- PCA with all the features
    print("\nSELECTED FEATURES\n")
    # Between the two conditions, regardless the duration of experiment (2 conditions, all durations)
    PCAOnDF(selected_df,
            target="condition",
            targets=["CHO", "DIO"],
            colors=["b", "r"],
            save_name="all_selected_features",
            save_in=save_in,
            three_D=True,
            )

    # # Between the two conditions, for each duration (2 conditions, 3 durations)
    # groupby_duration = selected_df.groupby("duration")
    # for duration, values in groupby_duration:
    #     # duration: str, values: pd_.DataFrame()
    #     PCAOnDF(values,
    #             target="condition",
    #             targets=["CHO", "DIO"],
    #             colors=["b", "r"],
    #             save_name=f"{duration}_selected_features",
    #             save_in=save_in,
    #             title=f" - {duration} Sample - selected features",
    #             plot_duration=False,
    #             three_D=True,
    #             )

    # -- K-means with all the features (2 conditions)

    # Between the two conditions, regardless the duration of experiment (2 conditions, all durations)
    kmeans = KmeansOnDF(selected_df,
                        target="condition",
                        nb_clusters=(2, 3, 4, 5),
                        elbow=False,
                        intracluster_var=False,
                        plot_bar=True,
                        save_name="all_selected_features",
                        save_in=save_in,
                        features_distribution=True,
                        )

    # # Between the two conditions, for each duration (2 conditions, 3 durations)
    # groupby_duration = selected_df.groupby("duration")
    # for duration, values in groupby_duration:
    #     kmeans = KmeansOnDF(values,
    #                         target="condition",
    #                         nb_clusters=(2,3,4,5),
    #                         elbow=False,
    #                         intracluster_var=True,
    #                         plot_bar=True,
    #                         save_name=f"{duration}_selected_features",
    #                         save_in=save_in,
    #                         title=f" - {duration} Sample - selected features",
    #                         duration=True,
    #                         )

    ## TODO: Random forests ?