Newer
Older
# 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 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,
Perform 2D or 3D PCA on the CHO-DIO dataframe.
Save explained variance ratio as csv, plot and save the PCAs.
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
# 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, 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, n_init=10, random_state=0)
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
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(nmpy.histogram(values['label'], bins=2))
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
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(nmpy.histogram(values['label'], bins=2))
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
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 = nmpy.arange(nb_cluster)
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
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 = nmpy.average([df.iloc[row, :] for row in soma_cluster], axis=0) # TODO .cluster_centers_
variance = sum([nmpy.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 = nmpy.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 = nmpy.triu(nmpy.ones_like(corr_matrix, dtype=nmpy.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 = nmpy.asarray(CHO[column])
DIO = cond_col_df.loc[cond_col_df['condition'] == "DIO"]
DIO = nmpy.asarray(DIO[column])
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
# 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_ = nmpy.asarray(CHO_[column])
DIO_ = nmpy.asarray(DIO_[column])
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
# 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(nmpy.triu(nmpy.ones(corr_matrix.shape), k=1).astype(nmpy.bool))
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
# 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 - nmpy.log(nmpy.sum(nmpy.sqrt(nmpy.multiply(Normalize(h1), Normalize(h2)))))
return h / nmpy.sum(h)
def 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 = nmpy.load("path.npy")
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
# 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,
# save_name="all_features_multiple_pop",
# save_in=save_in,
# features_distribution=False,
# )
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
# # 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,
)
## 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 RuntimeError("Only run the part until FeaturesStatistics included to generate df_drop_highly_corr_feat.csv, and then run the last part.")
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
## 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),
# intracluster_var=False,
# save_name="all_selected_features",
# save_in=save_in,
# )
# # 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 ?
if __name__ == "__main__":
#
Main()