<jupyter_start><jupyter_text>Moving Averages CINGULATE region  1) Imports<jupyter_code>import os
import glob
import sys
import json

p = os.path.abspath('../')
if p not in sys.path:
    sys.path.append(p)
print(sys.path)
import moving_averages as ma
import colorado as cld
import dico_toolbox as dtx
from tqdm import tqdm

import pandas as pd
import numpy as np
from scipy.spatial import distance

import plotly.graph_objects as go
import pickle
import matplotlib.pyplot as plt

from soma import aims

import torch

from sklearn.cluster import KMeans, SpectralClustering, AffinityPropagation
from sklearn import metrics
from sklearn.manifold import TSNE

import matplotlib.cm as cm

from hydra import initialize, initialize_config_module, initialize_config_dir, compose
from omegaconf import OmegaConf, DictConfig
def closest_distance(centroid, df):
    """Returns closest point to centroid of a given cluster
    /!\ central_point is part of the defined set of points
    Args:
        IN: centroid: tuple of coordinates of centroid
            df: pandas.DataFrame with index corresponding to subjects and as many columns as coordinates
        OUT: central_point: subject ID corresponding to closest point to the centroid
    """
    # df must have as many columns as dimensions to consider
    distances = {}
    # Get distances to centroid for each point
    for sub in list(df.index):
        pos = [df.loc[sub][k] for k in range(1, len(df.columns)+1)]
        distances[sub] = distance.euclidean(pos, centroid)
    # Get closest point to centroid
    central_point = min(distances, key=distances.get)
    return central_point
# Path the crops
path = '/neurospin/dico/data/deep_folding/current/crops/CINGULATE/mask/sulcus_based/2mm/centered_combined/tissier_2018/Rcrops/'

# Path to Tissier embeddings, created using loop_validate_clusterize_tissier.py
result_dir = '/host/volatile/jc225751/Runs/33_MIDL_2022_reviews/Output/t-0.1-analysis_tissier/n-004_o-4'

# Path to first half of hcp embeddings, created using loop_validate_clusterize.py
hcp_dir = '/host/volatile/jc225751/Runs/33_MIDL_2022_reviews/Output/t-0.1/n-004_o-4'

# Path To dir containing Tissier csv file with labels
input_dir = '/neurospin/dico/data/deep_folding/papers/midl2022/datasets/ACC_patterns'
# If true, meshes are saved as files
save_mesh = False

# if true, buckets are saved as files
save_bucket = False
def closest_distance(centroid, df):
    """Returns closest point to centroid of a given cluster
    /!\ central_point is part of the defined set of points
    Args:
        IN: centroid: tuple of coordinates of centroid
            df: pandas.DataFrame with index corresponding to subjects and as many columns as coordinates
        OUT: central_point: subject ID corresponding to closest point to the centroid
    """
    # df must have as many columns as dimensions to consider
    distances = {}
    # Get distances to centroid for each point
    for sub in list(df.index):
        pos = [df.loc[sub][k] for k in range(0, len(df.columns))]
        distances[sub] = distance.euclidean(pos, centroid)
    # Get closest point to centroid
    central_point = min(distances, key=distances.get)
    return central_point
def mscatter(x, y, ax=None, m=None, **kw):
    if not ax:
        ax = plt.gca()
    sc = ax.scatter(x, y, **kw)
    if (m is not None) and (len(m) == len(x)):
        paths = []
        for marker in m:
            if isinstance(marker, mmarkers.MarkerStyle):
                marker_obj = marker
            else:
                marker_obj = mmarkers.MarkerStyle(marker)
            path = marker_obj.get_path().transformed(
                marker_obj.get_transform())
            paths.append(path)
        sc.set_paths(paths)
    return sc<jupyter_output><empty_output><jupyter_text>Loading of subjects's embeddings and corresponding labelsDataFrame is obtained during analysis of beta-VAE latent space (2 dimensions) for CINGULATE crops, according to the following steps:- loading of trained model- encoding of test controls and asymmetry benchmark subjects<jupyter_code>if torch.cuda.is_available():
    device = "cuda:0"
print(device)
# We load the embeddings

val_embeddings_file = f"{result_dir}/train_val_embeddings.pt"
embeddings = torch.load(val_embeddings_file)
print(embeddings.shape)
# We load the labels, labels being subject names
with open(f"{result_dir}/train_val_filenames.json", 'r') as f:
    labels = json.load(f)
print(labels)
# We load the csv file containing labels (presence or not of parcacingulate sulcus)
tissier_labels = pd.read_csv(f"{input_dir}/tissier_labels.csv")
tissier_labels.head()
tissier_ACC = tissier_labels[['Sujet', 'ACC_right_2levels_new', 'ACC_left_2_levels_new']]
tissier_ACC.head()
tissier_ACC.ACC_right_2levels_new.replace(to_replace=dict(absent=0, present=1), inplace=True)
tissier_ACC.ACC_left_2_levels_new.replace(to_replace=dict(absent=0, present=1), inplace=True)
tissier_ACC.head()
kmeans = KMeans(n_clusters=3, random_state=0).fit(embeddings)
cluster_labels = kmeans.labels_
clusters_centroids = kmeans.cluster_centers_
print(f"clusters centroids = {clusters_centroids.shape}")
print(f"Average silhouette score: {metrics.silhouette_score(embeddings, cluster_labels)}")

dstrb_sub = pd.DataFrame(embeddings.numpy(), index=labels, columns=[1, 2, 3, 4])
dstrb_sub['cluster_lab'] = cluster_labels
dstrb_sub.head()
tissier_ACC.Sujet[0:3]
m2 = pd.merge(dstrb_sub, tissier_ACC, left_on=dstrb_sub.index, right_on=tissier_ACC.Sujet)
m2.head()
m2 = m2.set_index('Sujet')
m2 = m2.drop('key_0', 1)
m2.head()
dstrb_sub.shape
m2.shape<jupyter_output><empty_output><jupyter_text>Visualizations We first visualize on two first dimensions of the latent space<jupyter_code>arr = np.array([np.array([dstrb_sub[k][i] for k in m2.columns[0:2]]) for i in range(len(m2))])

color_dict = {0: 'blue', 1:'red'}
fig, ax = plt.subplots()

for g in np.unique([m2.ACC_right_2levels_new]):
    ix = np.where(m2.ACC_right_2levels_new == g)
    x = [arr[ix][k][0] for k in range(len(ix[0]))]
    y = [arr[ix][k][1] for k in range(len(ix[0]))]
    ax.scatter(x, y, c = color_dict[g], label = g)

plt.xlabel(f'dimension 1', fontsize=14)
plt.ylabel(f'dimension 2', fontsize=14)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>We then perform a t-SNE and visualize the results with labels<jupyter_code>m2_latent = m2.iloc[:,0:4]
m2_latent.head()
type(m2_latent)
embeddings = torch.tensor(m2_latent.values)
embeddings[:10]
tsne = TSNE(n_components=2, perplexity=5, init='pca', random_state=50)
X_tsne = tsne.fit_transform(embeddings)
lab_acc = list(m2['ACC_right_2levels_new'])
color_acc = [color_dict[lab] for lab in lab_acc]
lab_acc_left = list(m2['ACC_left_2_levels_new'])
color_acc_left = [color_dict[lab] for lab in lab_acc_left]
fig, ax = plt.subplots(1)
mscatter(X_tsne[:, 0], X_tsne[:, 1], c=color_acc_left, s=20, ax=ax)
plt.show()<jupyter_output><empty_output><jupyter_text>Get HCP data and put them together with Tissier data<jupyter_code># We load the hcp embeddings

val_embeddings_file = f"{hcp_dir}/train_val_embeddings.pt"
hcp_embeddings = torch.load(val_embeddings_file)
print(hcp_embeddings.shape)
# We load the labels, labels being subject names
with open(f"{hcp_dir}/train_val_filenames.json", 'r') as f:
    hcp_labels = json.load(f)
print(hcp_labels[:10])
tsne = TSNE(n_components=2, perplexity=30, init='pca', random_state=50)
X_tsne_hcp = tsne.fit_transform(hcp_embeddings)
fig, ax = plt.subplots(1)
mscatter(X_tsne_hcp[:, 0], X_tsne_hcp[:, 1], s=20, ax=ax)
plt.show()
hcp_df = pd.DataFrame(hcp_embeddings.numpy(), index=hcp_labels, columns=[1, 2, 3, 4])
hcp_df["cluster_lab"] = -1
hcp_df["ACC_right_2levels_new"] = -1
hcp_df["ACC_left_2_levels_new"] = -1
hcp_df.shape
hcp_df.index.names = ['Sujet']
hcp_df.shape
hcp_df.head()
m2.head()
all_df = hcp_df
all_df.head()
all_df.shape
all_df = all_df.append(m2)
all_df.shape
all_df.iloc[-1,:]
all_df_latent = all_df.iloc[:,0:4]
all_embeddings = torch.tensor(all_df_latent.values)
all_embeddings.shape
tsne = TSNE(n_components=2, perplexity=30, init='pca', random_state=50)
X_tsne_all = tsne.fit_transform(all_embeddings)
fig, ax = plt.subplots(1)
mscatter(X_tsne_all[:, 0], X_tsne_all[:, 1], s=20, ax=ax)
plt.show()
color_dict_all = {0: 'blue', 1:'red', -1: 'yellow'}
lab_all_acc = list(all_df['ACC_right_2levels_new'])
color_all_acc = [color_dict_all[lab] for lab in lab_all_acc]
lab_all_left_acc = list(all_df['ACC_left_2_levels_new'])
color_all_left_acc = [color_dict_all[lab] for lab in lab_all_left_acc]<jupyter_output><empty_output><jupyter_text>Correlated with right Acc paracingular labels<jupyter_code>fig, ax = plt.subplots(1)
mscatter(X_tsne_all[:, 0], X_tsne_all[:, 1], c=color_all_acc, s=20, ax=ax)
plt.show()<jupyter_output><empty_output><jupyter_text>Correlated with left Acc paracingular labels<jupyter_code>fig, ax = plt.subplots(1)
mscatter(X_tsne_all[:, 0], X_tsne_all[:, 1], c=color_all_left_acc, s=20, ax=ax)
plt.show()
all_df.head()
arr = np.array([np.array([all_df[k][i] for k in all_df.columns[0:2]]) for i in range(len(all_df))])

fig, ax = plt.subplots()

for g in np.unique([all_df.ACC_right_2levels_new]):
    ix = np.where(all_df.ACC_right_2levels_new == g)
    x = [arr[ix][k][0] for k in range(len(ix[0]))]
    y = [arr[ix][k][1] for k in range(len(ix[0]))]
    ax.scatter(x, y, c = color_dict_all[g], label = g)

plt.xlabel(f'dimension 1', fontsize=14)
plt.ylabel(f'dimension 2', fontsize=14)
plt.legend()
plt.show()
arr = np.array([np.array([all_df[k][i] for k in all_df.columns[[0,2]]]) for i in range(len(all_df))])

fig, ax = plt.subplots()

for g in np.unique([all_df.ACC_right_2levels_new]):
    ix = np.where(all_df.ACC_right_2levels_new == g)
    x = [arr[ix][k][0] for k in range(len(ix[0]))]
    y = [arr[ix][k][1] for k in range(len(ix[0]))]
    ax.scatter(x, y, c = color_dict_all[g], label = g)

plt.xlabel(f'dimension 1', fontsize=14)
plt.ylabel(f'dimension 3', fontsize=14)
plt.legend()
plt.show()
arr = np.array([np.array([all_df[k][i] for k in all_df.columns[1:3]]) for i in range(len(all_df))])

fig, ax = plt.subplots()

for g in np.unique([all_df.ACC_right_2levels_new]):
    ix = np.where(all_df.ACC_right_2levels_new == g)
    x = [arr[ix][k][0] for k in range(len(ix[0]))]
    y = [arr[ix][k][1] for k in range(len(ix[0]))]
    ax.scatter(x, y, c = color_dict_all[g], label = g)

plt.xlabel(f'dimension 2', fontsize=14)
plt.ylabel(f'dimension 3', fontsize=14)
plt.legend()
plt.show()
arr = np.array([np.array([all_df[k][i] for k in all_df.columns[[1,3]]]) for i in range(len(all_df))])

fig, ax = plt.subplots()

for g in np.unique([all_df.ACC_right_2levels_new]):
    ix = np.where(all_df.ACC_right_2levels_new == g)
    x = [arr[ix][k][0] for k in range(len(ix[0]))]
    y = [arr[ix][k][1] for k in range(len(ix[0]))]
    ax.scatter(x, y, c = color_dict_all[g], label = g)

plt.xlabel(f'dimension 2', fontsize=14)
plt.ylabel(f'dimension 4', fontsize=14)
plt.legend()
plt.show()
arr = np.array([np.array([all_df[k][i] for k in all_df.columns[2:4]]) for i in range(len(all_df))])

fig, ax = plt.subplots()

for g in np.unique([all_df.ACC_right_2levels_new]):
    ix = np.where(all_df.ACC_right_2levels_new == g)
    x = [arr[ix][k][0] for k in range(len(ix[0]))]
    y = [arr[ix][k][1] for k in range(len(ix[0]))]
    ax.scatter(x, y, c = color_dict_all[g], label = g)

plt.xlabel(f'dimension 3', fontsize=14)
plt.ylabel(f'dimension 4', fontsize=14)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Classification tests with Tissier database and cross-validation  SVM test<jupyter_code>from sklearn import preprocessing
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
from sklearn import svm
clf = make_pipeline(preprocessing.StandardScaler(), svm.SVC(C=1))
m2.head()
X = m2.iloc[:,0:4].values
y = m2.ACC_right_2levels_new.values
y
y_shuffle = m2.ACC_right_2levels_new.values.copy()
np.random.shuffle(y_shuffle)
y_shuffle
y
scores = cross_val_score(clf, X, y, cv=5)
print(scores)
print("%0.2f accuracy with a standard deviation of %0.2f" % (scores.mean(), scores.std()))
scores_shuffle = cross_val_score(clf, X, y_shuffle, cv=5)
print(scores_shuffle)
print("%0.2f accuracy with a standard deviation of %0.2f" % (scores_shuffle.mean(), scores_shuffle.std()))
# We do the scoring only on the best two dimensions of the latent
X_2_4 = X[:,[1,3]]
print(X[:3])
print(X_2_4[:3])
scores_2_4 = cross_val_score(clf, X_2_4, y, cv=5)
print(scores_2_4)
print("%0.2f accuracy with a standard deviation of %0.2f" % (scores_2_4.mean(), scores_2_4.std()))<jupyter_output>[0.55555556 0.875      0.57142857 0.71428571 0.71428571]
0.69 accuracy with a standard deviation of 0.12<jupyter_text>Logistic regression<jupyter_code>from sklearn.linear_model import LogisticRegression
clf_logit = LogisticRegression(solver='liblinear', random_state=0)
scores_logit = cross_val_score(clf_logit, X, y, cv=5)
print(scores_logit)
print("%0.2f accuracy with a standard deviation of %0.2f" % (scores_logit.mean(), scores_logit.std()))
scores_logit_shuffle = cross_val_score(clf_logit, X, y_shuffle, cv=5)
print(scores_logit_shuffle)
print("%0.2f accuracy with a standard deviation of %0.2f" % (scores_logit_shuffle.mean(), scores_logit_shuffle.std()))<jupyter_output>[0.55555556 0.625      0.57142857 0.42857143 0.57142857]
0.55 accuracy with a standard deviation of 0.07<jupyter_text>Gradient boosting classifier<jupyter_code>from sklearn.ensemble import GradientBoostingClassifier
clf_gdb = GradientBoostingClassifier(n_estimators=100, learning_rate=1.0, max_depth=1, random_state=0)
scores_gdb = cross_val_score(clf_gdb, X, y, cv=5)
print(scores_gdb)
print("%0.2f accuracy with a standard deviation of %0.2f" % (scores_gdb.mean(), scores_gdb.std()))<jupyter_output>[0.77777778 0.625      0.42857143 0.85714286 0.71428571]
0.68 accuracy with a standard deviation of 0.15
{'fraction_non_alphanumeric': 0.1316715753178651, 'fraction_numerical': 0.031822426549816825, 'mean_word_length': 5.38038496791934, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '327766', 'license': 'MIT', 'n_tokens_mistral': 5636, 'n_tokens_neox': 5328, 'n_words': 1218, 'path': 'contrastive/notebooks/10-analysis_tissier_dataset.ipynb', 'repo_name': 'neurospin-projects/2022_jchavas_cingulate_inhibitory_control', 'size': '13921'}