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_
import matplotlib.gridspec as gs_
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import scipy as si_
import seaborn as sb_

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 = ""
            ) -> list:
    '''
    Perform 2D PCA on the CHO-DIO dataframe.
    Print ratio variance and plot the PCA.
    '''
    # Separating the features from their conditions and durations
    all_target = pd_.DataFrame(df.loc[:, [target]].values, columns=[target])
    df_all = df.drop([target], 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)
    print(f"PCA explained variance ratio ({save_name}): ", pca.explained_variance_ratio_)

    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
    fig = pl_.figure(figsize=(8, 8))
    ax = fig.add_subplot(1, 1, 1)
    ax.set_xlabel('Principal Component 1', fontsize=15)
    ax.set_ylabel('Principal Component 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"D:\\MorganeNadal\\2_RESULTS\\Results_wo_erosion\\Features_analysis\\PCA_{save_name}.png")

    return pca.explained_variance_ratio_


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,
               title: str = "",
               ) -> 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.
    Returns kmeans.
    '''
    # Separating the features from their conditions and durations
    all_target = pd_.DataFrame(df.loc[:, [target]].values, columns=[target])
    df = df.drop([target], axis=1)

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

    # 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_.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)

        # Barplot
        if plot_bar:
            fig = pl_.figure(figsize=(8, 8))
            ax = fig.add_subplot(1, 1, 1)
            sb_.countplot(x="condition", hue="label", data=lab_cond_df, palette=sb_.color_palette("deep", n_colors=2))
            ax.set_title(f'Distribution of the clustering labels according to conditions{title}', fontsize=11)
            ax.grid()
            if save_name is not None:
                pl_.savefig(f"D:\\MorganeNadal\\2_RESULTS\\Results_wo_erosion\\Features_analysis\\Hist_Clustering_{save_name}.png")
                # pl_.show(block=True)
                # 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) -> list:
    '''
    Return the intracluster variance of a given cluster found by kmeans.
    '''
    var = []
    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)
        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 :", 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 = "",
                       ):
    '''
    Return the statistics allowing the user to choose the most relevant features to feed ML algorithms for ex.
    '''
    #
    # Overview of the basic stats on each columns of df
    description = df.describe()
    if save_in is not None:
        description.to_csv(f"{save_in}\\df_stat_description.csv")

    df_scalar = df.drop(["duration", "condition"], axis=1)
    df_cond = df.drop(["duration"], axis=1)
    df_groupby_cond = df.drop("duration", axis=1).groupby("condition")
    df_groupby_dur = df.groupby("duration")
    condition = pd_.DataFrame(df.loc[:, ["condition"]].values, columns=["condition"])

    # Overview of features distribution and correlation

    ## /!\ if too many features, error in seaborn display
    ## Even by splitting the data : too intense!

    # Plot heat map with correlation matrix btw features
    # Compute the correlation matrix
    print("Heat map with correlation matrix btw features")
    corr_matrix = df_scalar.corr().abs()
    # 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()

    # 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)
    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}\\df_drop_highly_corr_feat.csv")

    # Statistics for each features
    dict_ks = {}
    dict_wx = {}
    dict_ks_dur = {}
    dict_wx_dur = {}

    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()

        # Compare distribution between conditions (goodness of fit)
        print(f"Kolmogorov-Smirnov between conditions")
        ks = si_.stats.kstest(CHO[[column]], DIO[[column]])
        dict_ks[column] = ks

        # Compare median between conditions
        print(f"Wilcoxon signed-rank test between conditions")
        wx = si_.stats.wilcoxon(CHO[[column]], DIO[[column]])
        dict_wx[column] = wx

        # For each duration between conditions
        for duration, values in df_groupby_dur:
            dict_ks_dur[duration] = {}
            dict_wx_dur[duration] = {}

            duration_df = values.drop(["duration"])
            CHO_ = duration_df.loc[duration_df['condition'] == "CHO"]
            DIO_ = duration_df.loc[duration_df['condition'] == "DIO"]

            # Compare distribution
            print(f"Kolmogorov-Smirnov between conditions")
            ks2 = si_.stats.kstest(CHO[[column]], DIO[[column]])
            dict_ks_dur[duration][column] = ks2

            # Compare median
            print(f"Wilcoxon signed-rank test between conditions")
            wx2 = si_.stats.kstest(CHO[[column]], DIO[[column]])
            dict_wx_dur[duration][column] = wx2

        df_ks = pd_.DataFrame.from_dict()
        df_wx = pd_.DataFrame.from_dict()
        df_ks_dur = pd_.DataFrame.from_dict()
        df_wx_dur = pd_.DataFrame.from_dict()

        stat_tests_df = pd_.concatenate((df_ks, df_wx, df_ks_dur, df_wx_dur))


if __name__ == "__main__":
    #
    # 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)

    # path = sy_.argv[1]
    path = "D:\\MorganeNadal\\2_RESULTS\\Results_wo_erosion"

    df0 = pd_.read_csv(path,
                      # 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 (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
    # # Between the two conditions, regardless the duration of experiment (2 conditions, all durations)
    # df_all = df.drop(["duration"], axis=1)
    # all_pca = PCAOnDF(df_all, target="condition", targets=["CHO", "DIO"], colors=["b", "r"], save_name="all_features")
    #
    # # Between the two conditions, for each duration (2 conditions, 3 durations)
    # groupby_duration = df.groupby("duration")
    # duration_pca = []
    # for duration, values in groupby_duration:
    #     # print(duration, values.shape)
    #     # groupby_condition = values.groupby("condition")
    #     # for cond, val in groupby_condition:
    #     #     print(cond, val.shape)
    #     ## duration: str, values: pd_.DataFrame()
    #     duration_df = values.drop(["duration"], axis=1)
    #     pca = PCAOnDF(duration_df,
    #                   target="condition",
    #                   targets=["CHO", "DIO"],
    #                   colors=["b", "r"],
    #                   save_name=f"{duration}_features",
    #                   title=f" - {duration} Sample")
    #     duration_pca.append(pca)

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

    # # Between the two conditions, regardless the duration of experiment (2 conditions, all durations)
    # df_all = df.drop(["duration"], axis=1)
    # kmeans = KmeansOnDF(df_all,
    #                     target="condition",
    #                     nb_clusters=(2,),
    #                     elbow=False,
    #                     intracluster_var=True,
    #                     plot_bar=True,
    #                     save_name="all features"
    #                     )
    #
    # # Between the two conditions, for each duration (2 conditions, 3 durations)
    # groupby_duration = df.groupby("duration")
    # for duration, values in groupby_duration:
    #     duration_df = values.drop(["duration"], axis=1)
    #     kmeans = KmeansOnDF(duration_df,
    #                         target="condition",
    #                         nb_clusters=(2,),
    #                         elbow=False,
    #                         intracluster_var=True,
    #                         plot_bar=True,
    #                         save_name=f"{duration}_features",
    #                         title=f" - {duration} Sample",
    #                         )

    ## -- Various plots to analyse the data and find discriminant features by statistical analysis
    FeaturesStatistics(df,
                       save_in="D:\\MorganeNadal\\2_RESULTS\\Results_wo_erosion\\Features_analysis",
                       )
    ## TODO: Enter selected features here
    # selected_features = []
    # selected_df = df[selected_features]

    #
    # # -- PCA with selected features
    # # Between the two conditions, regardless the duration of experiment (2 conditions, all durations)
    # df_all = df.drop(["duration"], axis=1)
    # all_pca = PCAOnDF(df_all, target="condition", targets=["CHO", "DIO"], colors=["b", "r"], save_name="all_features")
    #
    # # Between the two conditions, for each duration (2 conditions, 3 durations)
    # groupby_duration = df.groupby("duration")
    # duration_pca = []
    # for duration, values in groupby_duration:
    #     # print(duration, values.shape)
    #     # groupby_condition = values.groupby("condition")
    #     # for cond, val in groupby_condition:
    #     #     print(cond, val.shape)
    #     ## duration: str, values: pd_.DataFrame()
    #     duration_df = values.drop(["duration"], axis=1)
    #     pca = PCAOnDF(duration_df,
    #                   target="condition",
    #                   targets=["CHO", "DIO"],
    #                   colors=["b", "r"],
    #                   save_name=f"{duration}_features",
    #                   title=f" - {duration} Sample")
    #     duration_pca.append(pca)

    # -- K-means with selected features (2 conditions)

    # # Between the two conditions, regardless the duration of experiment (2 conditions, all durations)
    # df_all = df.drop(["duration"], axis=1)
    # kmeans = KmeansOnDF(df_all,
    #                     target="condition",
    #                     nb_clusters=(2,),
    #                     elbow=False,
    #                     intracluster_var=True,
    #                     plot_bar=True,
    #                     save_name="all features"
    #                     )
    #
    # # Between the two conditions, for each duration (2 conditions, 3 durations)
    # groupby_duration = df.groupby("duration")
    # for duration, values in groupby_duration:
    #     duration_df = values.drop(["duration"], axis=1)
    #     kmeans = KmeansOnDF(duration_df,
    #                         target="condition",
    #                         nb_clusters=(2,),
    #                         elbow=False,
    #                         intracluster_var=True,
    #                         plot_bar=True,
    #                         save_name=f"{duration}_features",
    #                         title=f" - {duration} Sample",
    #                         )