utils_identifiability.py

"""
This file contains fast implementations for computing the Identifiability matrix of one
or two datasets, along with the PCA decomposition method proposed by E. Amico and J.Goni
in the paper "The quest for identifiability in human functional connectomes"

Authors: Andrea Santoro, edited by Ilaria Ricchi


# to integrate

def compute_MI_BOLD(i,j):
    global data_BOLD,mine
    mine.compute_score(data_BOLD[:,i], data_BOLD[:,j])
    mic=np.around(mine.mic(),3)
    return(i,j,mic)

"""

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from  multiprocessing import Pool
import matplotlib.pyplot as plt
import sys
from matplotlib.patches import Rectangle
import seaborn as sns
from matplotlib import colors
from matplotlib.gridspec import GridSpec
from scipy.io import loadmat
from scipy.spatial.distance import correlation as distcorr
from sklearn.feature_selection import mutual_info_regression as mi
from sklearn.metrics import mutual_info_score
from scipy.spatial.distance import cdist
from dtaidistance import dtw
from dtaidistance import dtw_ndim
from scipy import stats
import math
from math import sqrt
from scipy.linalg import sqrtm, logm, svd
import scipy.spatial.distance as dist
from utils_dtw import *
from minepy import cstats
from minepy import MINE

global mine

mine = MINE(alpha=0.6, c=15, est="mic_e")


## Style adjustements for plots
def plot_adjustments(ax=False, ticks_width=1.5, axis_width=2.5, labelsizereg=16):
    if ax == False:
        ax = plt.gca()
    #plt.rcParams['font.family'] = "PT Serif Caption"
    plt.rcParams['xtick.major.width'] = ticks_width
    plt.rcParams['ytick.major.width'] = ticks_width
    plt.rcParams['axes.linewidth'] = ticks_width
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["top"].set_linewidth(axis_width)
    ax.spines["right"].set_linewidth(axis_width)
    ax.spines["left"].set_linewidth(axis_width)
    ax.spines["bottom"].set_linewidth(axis_width)
    ax.tick_params(direction='out', length=4, width=axis_width + 0.5, colors='k',
                   grid_color='k', grid_alpha=0.5, axis='both')
    ax.tick_params(direction='out', which='minor', length=6, width=axis_width + 0.5, colors='k',
                   grid_color='k', grid_alpha=0.5, axis='both')
    ax.tick_params(direction='out', which='major', length=6, colors='k',
                   grid_color='k', grid_alpha=0.5, axis='both')
    ax.tick_params(axis='both', which='major', labelsize=labelsizereg)




#Taken from https://github.com/scikit-learn/scikit-learn/blob/f3f51f9b611bf873bd5836748647221480071a87/sklearn/utils/extmath.py
def svd_flip(u, v, u_based_decision=True):
    """Sign correction to ensure deterministic output from SVD.
    Adjusts the columns of u and the rows of v such that the loadings in the
    columns in u that are largest in absolute value are always positive.
    Parameters
    ----------
    u : ndarray
        u and v are the output of `linalg.svd` or
        :func:`~sklearn.utils.extmath.randomized_svd`, with matching inner
        dimensions so one can compute `np.dot(u * s, v)`.
    v : ndarray
        u and v are the output of `linalg.svd` or
        :func:`~sklearn.utils.extmath.randomized_svd`, with matching inner
        dimensions so one can compute `np.dot(u * s, v)`.
        The input v should really be called vt to be consistent with scipy's
        output.
    u_based_decision : bool, default=True
        If True, use the columns of u as the basis for sign flipping.
        Otherwise, use the rows of v. The choice of which variable to base the
        decision on is generally algorithm dependent.
    Returns
    -------
    u_adjusted, v_adjusted : arrays with the same dimensions as the input.
    """
    if u_based_decision:
        # columns of u, rows of v
        max_abs_cols = np.argmax(np.abs(u), axis=0)
        signs = np.sign(u[max_abs_cols, range(u.shape[1])])
        u *= signs
        v *= signs[:, np.newaxis]
    else:
        # rows of v, columns of u
        max_abs_rows = np.argmax(np.abs(v), axis=1)
        signs = np.sign(v[range(v.shape[0]), max_abs_rows])
        u *= signs
        v *= signs[:, np.newaxis]
    return u, v


#https://stackoverflow.com/questions/71844846/is-there-a-faster-way-to-get-correlation-coefficents
def pairwise_correlation(A, B):
    '''
    Compute the Pearson correlation coefficient between two matrices

    input: A -> numpy 2D array
           B -> numpy 2D array

    output: matrix encoding all the Pearson correlation'''

    
    am = A - np.mean(A, axis=0, keepdims=True)
    bm = B - np.mean(B, axis=0, keepdims=True)
    return am.T @ bm /  (np.sqrt(
        np.sum(am**2, axis=0,
               keepdims=True)).T * np.sqrt(
        np.sum(bm**2, axis=0, keepdims=True)))

# def riemannian_distance(fc1, fc2):
#     """
#     Compute the Riemannian distance between two FC matrices, ensuring diagonal values are ones.
#     """
#     fc1 = fc1 @ fc1.T  # Making it symmetric positive-definite
#     np.fill_diagonal(fc1, 1)  # Ensuring diagonal is ones
#     fc2 = fc2 @ fc2.T  # Making it symmetric positive-definite
#     np.fill_diagonal(fc1, 1)  # Ensuring diagonal is ones
    
#     # np.fill_diagonal(fc1, 1)
#     # np.fill_diagonal(fc2, 1)
    
#     sqrt_fc1 = sqrtm(fc1)
#     sqrt_fc2 = sqrtm(fc2)
    
#     log_diff = logm(sqrt_fc1 @ np.linalg.inv(fc2) @ sqrt_fc1)
#     return np.linalg.norm(log_diff, 'fro')

def riemannian_distance(P0, P1):
    """
    Compute the Riemannian distance between two symmetric positive-definite FC matrices.
    
    Parameters:
    - P0: ndarray (n x n), Reference FC matrix
    - P1: ndarray (n x n), Target FC matrix

    Returns:
    - distance: float, Riemannian distance
    """
    # Compute the matrix square root inverse of P0
    P0_sqrt_inv = np.linalg.inv(sqrtm(P0))  # Equivalent to (P0)^(-1/2)

    # Compute the congruence transform
    M = P0_sqrt_inv @ P1 @ P0_sqrt_inv  # (P0)^(-1/2) * P1 * (P0)^(-1/2)

    # Compute the matrix logarithm using SVD
    U, S, Vh = svd(M)  # Singular Value Decomposition
    S_log = np.diag(np.log(S))  # Log of singular values
    logM = U @ S_log @ Vh  # Reconstruct the log matrix

    # Compute the Frobenius norm
    distance = np.linalg.norm(logM, 'fro')

    return distance

def riemannian_similarity(fc1, fc2, alpha=1.0):
    """
    Compute the equivalent of correlation using the exponential of the negative Riemannian distance.
    """
    d_r = riemannian_distance(fc1, fc2)
    return 1-d_r/np.max(d_r) #np.exp(-alpha * d_r)

def compute_riemannian_Imat(fc_list1, fc_list2):
    """
    Compute a similarity matrix for lists of FC matrices.
    """
    n = len(fc_list1)
    m = len(fc_list2)
    assert n == m

    # FC_list1 = [fc @ fc.T for fc in fc_list1]  # Making them symmetric positive-definite
    # FC_list2 = [fc @ fc.T for fc in fc_list2]  # Making them symmetric positive-definite
    
    Imat = np.zeros((n, n))
    print(fc_list1)
    print(fc_list1[0].shape)

    print(fc_list2)
    print(fc_list2[0].shape)
    for i in range(n):
        for j in range(m):
            Imat[i, j] = riemannian_similarity(fc_list1[i], fc_list2[j])
    
    return Imat
    

def cohen_d(x,y):
    return (np.mean(x) - np.mean(y)) / sqrt((np.std(x, ddof=1) ** 2 + np.std(y, ddof=1) ** 2) / 2.0)

def compute_Iscores(Imat):
    '''
    Function to compute the Idiff, Iself and Iothers from an identifiability matrix Imat

    input: Identiability matrix (nsubjs x nsubjs)

    output: Differential Identifiability score (%),
            Iself
            Iothers
    '''
    Iself=np.mean(Imat.diagonal())
    Iothers=np.mean(Imat[~np.eye(Imat.shape[0],dtype=bool)])
    Idiff=(Iself-Iothers)*100

    return(Idiff,Iself,Iothers)

def from_triangle2matrix(upper_values, n):
    """Reconstruct a symmetric n x n matrix from its upper triangle values."""
    full_matrix = np.zeros((n, n))
    u, v = np.triu_indices(n, k=1)  # Get upper triangle indicaes
    full_matrix[u, v] = upper_values  # Fill upper triangle
    full_matrix[v, u] = upper_values  # Copy to lower triangle
    return full_matrix

def computing_FCs_halved(data,trimming_length=0,metric="pairwise_correlation",inter_bet_con=None):
    ''' 
    Function for computing Test and Retest  Functional Connectomes from a 
    dictionary in the form {subjID: fMRI data} by splitting each fMRI
    recording in two halves 

    input: data -> dictionary in the form {subjID: fMRI data}
           trimming_length -> default 0 (no trimming), if different from 0 then trim
                              all the recordings to the same length
           metric -> pairwise correlation is the default, it accepts: MI, TE, distcorr

    output: a 3D array of size (nsubjects x 2 x nedges) reporting the 
            flattened upper triangular matrix of functional connectomes 
            constructed when splitting each fMRI recording in two halves
    '''

    nsubjects=len(data)
    n_ROIs=np.shape(list(data.values())[0])[0]
    
    u,v=np.triu_indices(n=n_ROIs,k=1,m=n_ROIs)
    FC_list=[];
    
    ##Computing the functional connectomes splitting the TC in half
    for x in data.values():
        x=np.array(x)
        if trimming_length!=0:
            TC=x[:,:trimming_length]    
        else:
            TC=x
        n_timepoints=np.shape(TC)[1]
        half_T=int(n_timepoints/2)
        
        if metric == "pairwise_correlation":
            FC_1=np.nan_to_num(pairwise_correlation(TC[:,:half_T].T,TC[:,:half_T].T))
            FC_2=np.nan_to_num(pairwise_correlation(TC[:,half_T:].T,TC[:,half_T:].T))
            
        elif metric == "MI":
            FC_1 = np.zeros((n_ROIs,n_ROIs))
            FC_2 = np.zeros((n_ROIs,n_ROIs))
            
            for i in range(n_ROIs):
                FC_1[i,:] = mi(TC[:,:half_T].T,TC[i,:half_T], n_neighbors=7)
                FC_2[i,:] = mi(TC[:,half_T:].T,TC[i,half_T:], n_neighbors=7)

        elif metric == "dist_correlation":
            FC_1 = np.zeros((n_ROIs,n_ROIs))
            FC_2 = np.zeros((n_ROIs,n_ROIs))
            for i in range(n_ROIs):
                for j in range(n_ROIs):
                    FC_1[i,j] = distcorr(TC[i,:half_T].T,TC[j,:half_T].T)
                    FC_2[i,j] = distcorr(TC[i,half_T:].T,TC[j,half_T:].T)  
        elif metric == "dtw":
            dist_1 = dtw.distance_matrix_fast(TC[:,:half_T],window=2)
            dist_2 = dtw.distance_matrix_fast(TC[:,half_T:],window=2)
            dist_1_n = dist_1/np.max(dist_1)
            dist_2_n = dist_2/np.max(dist_2)
            FC_1 = 1-dist_1_n
            FC_2 = 1-dist_2_n
            
        FC_list.append([FC_1[u,v],FC_2[u,v]])
        
    FC_list=np.array(FC_list)
    return FC_list

def computing_FCs_halved_concat(data1,data2,data,trimming_length=0):
    ''' 
    Function for computing Test and Retest  Functional Connectomes from a 
    dictionary in the form {subjID: fMRI data} by splitting each fMRI
    recording in two halves 

    
    input: data1/2 -> dictionary in the form {subjID: fMRI data}
           trimming_length -> default 0 (no trimming), if different from 0 then trim
                              all the recordings to the same length

    output: a 3D array of size (nsubjects x 2 x nedges) reporting the 
            flattened upper triangular matrix of functional connectomes 
            constructed when splitting each fMRI recording in two halves
    '''
    assert len(data1)== len(data2)
    
    
    nsubjects=len(data1)
    n_ROIs1=np.shape(list(data1.values())[0])[0]
    n_ROIs2=np.shape(list(data2.values())[0])[0]
    u1,v1=np.triu_indices(n=n_ROIs1,k=1,m=n_ROIs1)
    u2,v2=np.triu_indices(n=n_ROIs2,k=1,m=n_ROIs2)
    
    u = np.concatenate([u1,u2+n_ROIs1])
    v = np.concatenate([v1,v2+n_ROIs1])
    FC_list=[];
    
    ##Computing the functional connectomes splitting the TC in half
    for x in data.values():
        x=np.array(x)
        if trimming_length!=0:
            TC=x[:,:trimming_length]    
        else:
            TC=x
        n_timepoints=np.shape(TC)[1]
        half_T=int(n_timepoints/2)
        
        FC_1=np.nan_to_num(pairwise_correlation(TC[:,:half_T].T,TC[:,:half_T].T))
        FC_2=np.nan_to_num(pairwise_correlation(TC[:,half_T:].T,TC[:,half_T:].T))
        FC_list.append([FC_1[u,v],FC_2[u,v]])
    FC_list=np.array(FC_list)
    return FC_list


def computing_FCs(data,trimming_length=0,metric="pairwise_correlation",inter_bet_con=None,W=1):
    '''
    Function for computing Functional Connectomes/Correlation Matrices from a 
    dictionary in the form {subjID: fMRI data}

    input: data -> dictionary in the form {subjID: fMRI data}
           trimming_length -> default 0 (no trimming), if different from 0 then trim
                              all the recordings to the same length

    output: a 2D array of size (nsubjects x nedges) reporting the 
            flattened upper triangular matrix of functional connectomes 
            constructed from the fMRI recordings of each subject'''

    nsubjects=len(data)
    n_ROIs=np.shape(list(data.values())[0])[0]
    
    u,v=np.triu_indices(n=n_ROIs,k=1,m=n_ROIs)
    FC_list=[];
    
    ##Computing the functional connectomes for each subject
    for x in data.values():
        x=np.array(x)
        FC = np.zeros((n_ROIs,n_ROIs))
        
        if trimming_length!=0:
            TC=x[:,:trimming_length]    
        else:
            TC=x
            
        if metric == "pairwise_correlation":
            FC=np.nan_to_num(W*pairwise_correlation(TC.T,TC.T))
        elif metric == "MI":
            for i in range(n_ROIs):
                FC[i,:] = mi(TC.T,TC[i,:], n_neighbors=7)
        elif metric == "dist_correlation":
            for i in range(n_ROIs):
                for j in range(n_ROIs):
                    FC[i,j] = distcorr(TC[i,:].T,TC[j,:].T)
        elif metric == "dtw":
            dist = dtw.distance_matrix_fast(TC)
            dist_n = dist/np.max(dist)
            FC = np.exp( - dist_n)

        if inter_bet_con is not None:
            br = inter_bet_con['br']
            sc = inter_bet_con['sp']
            FC_list.append(FC[br:br+sc,:br].flatten())
        else:
            FC_list.append(FC[u,v])
            
    FC_list=np.array(FC_list)
    return (FC_list)


def compute_Imat_from_single_session_TC(data,trimming_length=0,nodes=None,comp_metric='pearson_correlation'):
    '''
    Function to compute Imat, Idiff, Iself and Iothers from a dictionary in the form
    {subjID: fMRI data}.

    input: data -> dictionary computing_FCs_halved in the form {subjID: fMRI data}
           trimming_length -> default 0 (no trimming), if different from 0 then trim
                              all the recordings to the same length

    output: Identifiability matrix,
            Differential Identifiability score (%),
            Iself
            Iothers
    '''
    nsubjects=len(data)
    n_ROIs=np.shape(list(data.values())[0])[0]
    u,v=np.triu_indices(n=n_ROIs,k=1,m=n_ROIs)
    FC_list=computing_FCs_halved(data,trimming_length)
    ##Compute the Identifiability matrix
    if comp_metric == 'pearson_correlation':
        I_mat=pairwise_correlation(FC_list[:,0].T,FC_list[:,1].T)
    elif comp_metric == 'riemannian_sim':
        FC_list_mat1 = [from_triangle2matrix(fc,n_ROIs)  for fc in FC_list[:,0]]
        FC_list_mat2 = [from_triangle2matrix(fc,n_ROIs)  for fc in FC_list[:,1]]
        I_mat=compute_riemannian_Imat(FC_list_mat1,FC_list_mat2)
    elif comp_metric == 'MI':
        I_mat = np.zeros((nsubjects,nsubjects))
        for i in range(nsubjects):
            for j in range(nsubjects):
                I_mat[i,j] = mutual_info_score(FC_list[:,0][i,:],FC_list[:,0][j,:])
        
    Idiff,Iself,Iothers=compute_Iscores(I_mat)
    return(I_mat,Idiff,Iself,Iothers)


def compute_Imat_from_TC_twodatasets(data1,data2,trimming_length=0,comp_metric='pearson_correlation'):
    '''
    Function to compute Imat, Idiff, Iself and Iothers from two dictionaries data1 and 
    data2 in the form {subjID: fMRI data}, data1 is Test, data2 is Retest 

    input: data1 -> dictionary in the form {subjID: fMRI data}
           data2 -> dictionary in the form {subjID: fMRI data}
           trimming_length -> default 0 (no trimming), if different from 0 then trim
                              all the recordings to the same length

    output: Identifiability matrix,
            Differential Identifiability score (%),
            Iself
            Iothers
    '''
    nsubjects=len(data1)
    n_ROIs=np.shape(list(data1.values())[0])[0]
    u,v=np.triu_indices(n=n_ROIs,k=1,m=n_ROIs)
    FC1_list=computing_FCs(data1,trimming_length)
    FC2_list=computing_FCs(data2,trimming_length)
    ##Analogous command of np.array(list(zip(FC1_list,FC2_list)))) but slightly faster
    FC_list=np.stack((FC1_list,FC2_list), axis=1)
    ##Compute the Identifiability matrix
    if comp_metric == 'pearson_correlation':
        I_mat=pairwise_correlation(FC_list[:,0].T,FC_list[:,1].T)
    elif comp_metric == 'riemannian_sim':
        FC_list_mat1 = [from_triangle2matrix(fc,n_ROIs)  for fc in FC_list[:,0]]
        FC_list_mat2 = [from_triangle2matrix(fc,n_ROIs)  for fc in FC_list[:,1]]
        I_mat=compute_riemannian_Imat(FC_list_mat1,FC_list_mat2)
    elif comp_metric == 'MI':
        I_mat = np.zeros((nsubjects,nsubjects))
        for i in range(nsubjects):
            for j in range(nsubjects):
                I_mat[i,j] = mutual_info_score(FC_list[:,0][i,:],FC_list[:,0][j,:])
        
    Idiff,Iself,Iothers=compute_Iscores(I_mat)
    return(I_mat,Idiff,Iself,Iothers)    

def retrieve_n(num_edges):
    """Retrieve n from the number of upper triangular edges."""
    n = (1 + math.sqrt(1 + 8 * num_edges)) / 2
    return int(n) if n.is_integer() else None  # Ensure n is a whole number


def compute_Imat_from_single_session_FCs(FC_list,comp_metric='pearson_correlation'):
    '''
    Function to compute Imat, Idiff, Iself and Iothers from a 3D array of FCs
    of the form (nsubjs x 2 x nedges)

    input: data -> numpy array with a shape: nsubjs x 2 x nedges(upper triangular FC)
    
    output: Identifiability matrix,
            Differential Identifiability score (%),
            Iself
            Iothers
    '''
    nsubjects=len(FC_list)
    FC_list=np.array(FC_list)
    num_edges = len(FC_list[0][0])
    n_ROIs = retrieve_n(num_edges)
    ##Computing the Identifiability matrix
    if comp_metric == 'pearson_correlation':
        I_mat=pairwise_correlation(FC_list[:,0].T,FC_list[:,1].T)   
    elif comp_metric == 'riemannian_sim':
        FC_list_mat1 = [from_triangle2matrix(fc,n_ROIs)  for fc in FC_list[:,0]]
        FC_list_mat2 = [from_triangle2matrix(fc,n_ROIs)  for fc in FC_list[:,1]]
        I_mat=compute_riemannian_Imat(FC_list_mat1,FC_list_mat2)
    elif comp_metric == 'MI':
        I_mat = np.zeros((nsubjects,nsubjects))
        for i in range(nsubjects):
            for j in range(nsubjects):
                I_mat[i,j] = mutual_info_score(FC_list[:,0][i,:],FC_list[:,0][j,:])
        
    Idiff,Iself,Iothers=compute_Iscores(I_mat)
    return(I_mat,Idiff,Iself,Iothers)


def compute_Imat_from_FCs_twodatasets(FC1_list,FC2_list, comp_metric='pearson_correlation'):
    '''
    Function to compute Imat, Idiff, Iself and Iothers from two 2D arrays of size (nsubjects x nedges)
    reporting the flattened upper triangular matrix of functional connectomes 
    constructed from the fMRI recordings of each subject

    input: FC1_list -> a 2D array of size (nsubjects x nedges) with flattened FCs for TEST
           FC2_list -> a 2D array of size (nsubjects x nedges) with flattened FCs for RETEST

    output: Identifiability matrix,
            Differential Identifiability score (%),
            Iself
            Iothers
    '''
    ##Analogous command of np.array(list(zip(FC1_list,FC2_list)))) but slightly faster
    nsubjects=FC1_list.shape[0]
    FC_list=np.stack((FC1_list,FC2_list), axis=1)
    num_edges = len(FC_list[0][0])
    n_ROIs = retrieve_n(num_edges)
    # weighted (for KCL dataset)
    if comp_metric == 'pearson_correlation':
        I_mat=pairwise_correlation(FC_list[:,0].T,FC_list[:,1].T)
    elif comp_metric == 'riemannian_sim':
        FC_list_mat1 = [from_triangle2matrix(fc,n_ROIs)  for fc in FC1_list]
        FC_list_mat2 = [from_triangle2matrix(fc,n_ROIs)  for fc in FC2_list]
        I_mat=compute_riemannian_Imat(FC_list_mat1,FC_list_mat2)
    elif comp_metric == 'MI':
        I_mat = np.zeros((nsubjects,nsubjects))
        for i in range(nsubjects):
            for j in range(nsubjects):
                I_mat[i,j] = mutual_info_score(FC1_list[i,:],FC1_list[j,:])
        
    Idiff,Iself,Iothers=compute_Iscores(I_mat)
    
    return(I_mat,Idiff,Iself,Iothers) 



def compute_Imat_from_listFCs(data_FC,comp_metric='pearson_correlation'):
    '''
    Function to compute Imat, Idiff, Iself and Iothers from a 2D array of FCs
    of the form (2*nsubjs x nedges), alternating Test, Retest as odd and even
    rows 

    input: data -> numpy array with a shape: 2*nsubjs x nedges(upper triangular FC)

    output: Identifiability matrix,
            Differential Identifiability score (%),
            Iself
            Iothers
    '''
    nsubj_TR,nlinks=np.shape(data_FC)
    ##nsubj_TR is double of the number of subjects since it contains test and retes
    nsubj=int(nsubj_TR/2)##These are all the number of subjects
    #Reshape the matrix as a 3D array
    data_FC=data_FC.reshape((nsubj,2,nlinks))
    n_ROIs = retrieve_n(nlinks)
    if comp_metric == 'pearson_correlation':
        Imat=pairwise_correlation(data_FC[:,0].T,data_FC[:,1].T)
    elif comp_metric == 'riemannian_sim':
        FC_list_mat1 = [from_triangle2matrix(fc,n_ROIs)  for fc in data_FC[:,0]]
        FC_list_mat2 = [from_triangle2matrix(fc,n_ROIs)  for fc in data_FC[:,1]]
        I_mat=compute_riemannian_Imat(FC_list_mat1,FC_list_mat2)
    elif comp_metric == 'MI':
        I_mat = np.zeros((nsubj,nsubj))
        for i in range(nsubj):
            for j in range(nsubj):
                I_mat[i,j] = mutual_info_score(FC_list[:,0][i,:],FC_list[:,0][j,:])
      

    Idiff,Iself,Iothers=compute_Iscores(Imat)
    return(Imat,Idiff,Iself,Iothers)


def fit_kcomponents(U,S,Vt,all_connectomes,mean_,k_components):
    '''
    Function that reconstruct the matrix "all_connectomes"
    starting from the first k_components of the pca 
    (and decomposition U,S,Vt) passed as input 

    input: U,S,Vt-> Decomposition obtained from the command svd decomposition of all connectomes
           all_connectomes -> 2d numpy array containing all the flattened FC (upper triangular)
                             that was passed as input for the pca._fit(all_connectomes)
           mean_ -> mean of all the connectomes across rows (i.e. np.mean(all_connectomes,axis=0))
           k_components -> integer from 1 to Ncomponents max 

    output: 2D array containing the reconstructed numpy array of "all_connectomes" using only
           k_components
    ''' 
    
    UU = U[:, : k_components]
    UU =UU* S[: k_components]
    reconstructed=np.dot(UU, Vt[:k_components]) + mean_
    return(reconstructed)

def fit_k_sel_components(U,S,Vt,all_connectomes,mean_,sel_components):
    '''
    Function that reconstruct the matrix "all_connectomes"
    starting from the first k_components of the pca 
    (and decomposition U,S,Vt) passed as input 

    input: U,S,Vt-> Decomposition obtained from the command svd decomposition of all connectomes
           all_connectomes -> 2d numpy array containing all the flattened FC (upper triangular)
                             that was passed as input for the pca._fit(all_connectomes)
           mean_ -> mean of all the connectomes across rows (i.e. np.mean(all_connectomes,axis=0))
           sel_components -> integer indexes of the components to consider

    output: 2D array containing the reconstructed numpy array of "all_connectomes" using only
           k_components
    ''' 
    
    UU = U[:, sel_components]
    UU =UU* S[sel_components]
    reconstructed=np.dot(UU, Vt[sel_components]) + mean_
    return(reconstructed)


def compute_Idiff_PCA_one_dataset_from_TC(data,trimming_length=0,comp_metric="pearson_correlation"):
    '''
    Compute the Idiff values as a function of the number of components
    as proposed in the paper by Amico and Goni "The quest for identifiability
    in human functional connectomes"

    input: data -> dictionary in the form {subjID: fMRI data}
           trimming_length -> default 0 (no trimming), if different from 0 then trim
                              all the recordings to the same length

    output: numpy array containing the following rows:
            (k_components,Differential Identifiability score (%))
    '''

    nsubjects=len(data)
    list_FCs=computing_FCs_halved(data,trimming_length=trimming_length)
    all_connectomes=list_FCs.reshape((nsubjects*2,len(list_FCs[0][0])))
    N_total=len(all_connectomes)
    mean_ = np.mean(all_connectomes, axis=0)
    all_connectomes -=mean_
    U, S, Vt= np.linalg.svd(all_connectomes, full_matrices=False)
    U, Vt = svd_flip(U, Vt)
    # pca = PCA(n_components=N_total)
    # U, S, Vt=pca._fit(all_connectomes)

    ## Doing the PCA reconstruction keeping k components
    list_Idiff=np.array([(k_components,
      compute_Imat_from_listFCs(fit_kcomponents(U,S,Vt,all_connectomes,mean_,k_components),comp_metric)[1])
      for k_components in range(1,len(all_connectomes)+1)])

#    #The previous lines are equivalent to these
#     list_Idiff=[]    
#     for k_components in range(1,len(all_connectomes)+1):
#         reconstructed_connectomes=fit_kcomponents(pca,U,S,Vt,all_connectomes,k_components)
#         _,current_Idiff,_,_=compute_Imat_from_listFCs(reconstructed_connectomes)
#         list_Idiff.append([k_components,current_Idiff])
#     list_Idiff=np.array(list_Idiff)

    return(list_Idiff)


def compute_Idiff_PCA_one_dataset_from_FCs(list_FCs,optimal=False,comp_metric="pearson_correlation"):
    '''
    Compute the Idiff values as a function of the number of components
    as proposed in the paper by Amico and Goni "The quest for identifiability
    in human functional connectomes"

    input: list_FCs -> a 3D array of size (nsubjects x 2 x nedges) reporting the 
                       flattened upper triangular matrix of functional connectomes 
                       constructed when splitting each fMRI recording in two halves
           optimal -> bool value, if optimal is set to True, the function returns
                      the optimial Identifiability matrix obtained from the PCA
                      decomposition method
    output: numpy array containing the following rows:
            (k_components,Differential Identifiability score (%))
            if optimal is True, the returns also the optimal Identifiability matrix
    '''

    nsubjects=len(list_FCs)
    all_connectomes=list_FCs.reshape((nsubjects*2,len(list_FCs[0][0])))
    N_total=len(all_connectomes)
    mean_ = np.mean(all_connectomes, axis=0)
    all_connectomes -=mean_
    U, S, Vt= np.linalg.svd(all_connectomes, full_matrices=False)
    U, Vt = svd_flip(U, Vt)
    # pca = PCA(n_components=N_total)
    # U, S, Vt=pca._fit(all_connectomes)

    ## Doing the PCA reconstruction keeping k components
    list_Idiff=np.array([(k_components,
      compute_Imat_from_listFCs(fit_kcomponents(U,S,Vt,all_connectomes,mean_,k_components),comp_metric)[1])
      for k_components in range(1,len(all_connectomes)+1)])

    
    # opt_Imat = None
    # if optimal:
        
    #     print(list_Idiff)
    #     ncomp_opt=int(np.argmax(list_Idiff[:,1]))
    #     ncomp_opt=int(list_Idiff[ncomp_opt,0])
    #     opt_Imat,_,_,_=compute_Imat_from_listFCs(fit_kcomponents(U,S,Vt,all_connectomes,mean_,ncomp_opt))


    opt_Imat = None
    ncomp_opt= None    
    print(optimal, type(optimal))
    if isinstance(optimal,int):
        ncomp_opt=optimal
        
    elif optimal:
        print(list_Idiff)
        ncomp_opt=int(np.argmax(list_Idiff[:,1]))
        ncomp_opt=int(list_Idiff[ncomp_opt,0])
    
    if optimal != None:
        opt_Imat,_,_,_=compute_Imat_from_listFCs(fit_kcomponents(U,S,Vt,all_connectomes,mean_,ncomp_opt),comp_metric)
    

    return(list_Idiff,opt_Imat)


def compute_Idiff_PCA_two_datasets_from_TC(data1,data2,trimming_length=0,comp_metric="pearson_correlation"):
    '''
    Compute the Idiff values as a function of the number of components
    as proposed in the paper by Amico and Goni "The quest for identifiability
    in human functional connectomes"

    input: data1 -> dictionary in the form {subjID: fMRI data} %% TEST
           data2 -> dictionary in the form {subjID: fMRI data} %% RETEST
           trimming_length -> default 0 (no trimming), if different from 0 then trim
                              all the recordings to the same length

    output: numpy array containing the following rows:
            (k_components,Differential Identifiability score (%))

    '''

    nsubjects=len(data1)
    n_ROIs=np.shape(list(data1.values())[0])[0]
    u,v=np.triu_indices(n=n_ROIs,k=1,m=n_ROIs)
    FC1_list=computing_FCs(data1,trimming_length)
    FC2_list=computing_FCs(data2,trimming_length)
    ##Analogous command of np.array(list(zip(FC1_list,FC2_list)))) but slightly faster
    list_FCs=np.stack((FC1_list,FC2_list), axis=1)
    #list_FCs=computing_FCs_halved(data,trimming_length=trimming_length)
    all_connectomes=list_FCs.reshape((nsubjects*2,len(list_FCs[0][0])))
    N_total=len(all_connectomes)
    # pca = PCA(n_components=N_total)
    # U, S, Vt=pca._fit(all_connectomes)
    mean_ = np.mean(all_connectomes, axis=0)
    all_connectomes -=mean_
    U, S, Vt= np.linalg.svd(all_connectomes, full_matrices=False)
    U, Vt = svd_flip(U, Vt)

    ## Doing the PCA reconstruction keeping k components
    list_Idiff=np.array([(k_components,
      compute_Imat_from_listFCs(fit_kcomponents(U,S,Vt,all_connectomes,mean_,k_components),comp_metric)[1])
      for k_components in range(1,len(all_connectomes)+1)])

    return(list_Idiff)


def Idiff_from_FC_parallel(k_components,comp_metric="pearson_correlation"):
    reconstructed_connectome=fit_kcomponents(U,S,Vt,list_connectomes,mean_,k_components)
    _,Idiff,_,_=compute_Imat_from_listFCs(reconstructed_connectome,comp_metric)
    return(k_components,Idiff)

def create_pca_approach_parallel(all_connectomes,N_total):
    global list_connectomes, U, S, Vt,mean_
    list_connectomes=all_connectomes
    mean_ = np.mean(list_connectomes, axis=0)
    list_connectomes -= mean_
    U, S, Vt= np.linalg.svd(list_connectomes, full_matrices=False)
    U, Vt = svd_flip(U, Vt)

def compute_Idiff_PCA_two_datasets_from_FCs(FC1_list,FC2_list,parallel=None, optimal=False,comp_metric="pearson_correlation"):
    '''
    Compute the Idiff values as a function of the number of components
    as proposed in the paper by Amico and Goni "The quest for identifiability
    in human functional connectomes" from two lists of FCs

    input: FC1_list -> a 2D array of size (nsubjects x nedges) with flattened FCs for TEST
           FC2_list -> a 2D array of size (nsubjects x nedges) with flattened FCs for RETEST

    output: numpy array containing the following rows:
            (k_components,Differential Identifiability score (%))
            
    '''

    nsubjects=len(FC1_list)
    ##Analogous command of np.array(list(zip(FC1_list,FC2_list)))) but slightly faster
    list_FCs=np.stack((FC1_list,FC2_list), axis=1)
    all_connectomes=list_FCs.reshape((nsubjects*2,len(FC1_list[0])))
    N_total=len(all_connectomes)
    varexp = None
    if parallel==None:
        #Doing the PCA
        mean_ = np.mean(all_connectomes, axis=0)
        all_connectomes -=mean_
        U, S, Vt= np.linalg.svd(all_connectomes, full_matrices=False)
        U, Vt = svd_flip(U, Vt)
        ## Doing the PCA reconstruction keeping k components
        ## the original here
        list_Idiff=np.array([(k_components,
          compute_Imat_from_listFCs(fit_kcomponents(U,S,Vt,all_connectomes,mean_,k_components),comp_metric)[1])
          for k_components in range(1,len(all_connectomes)+1)])
        
        varexp = (np.diag(S)**2 / np.sum(np.diag(S)**2))
    else:
        ###Parallel processing is slower than sequential (due to to )
        results=[]
        pool = Pool(processes=parallel, initializer=create_pca_approach_parallel,
                initargs=(all_connectomes, N_total))
        
        for k_components in range(1,len(all_connectomes)+1):
            results.append(pool.apply_async(Idiff_from_FC_parallel,args=(k_components,)))
        pool.close()
        pool.join()
        list_Idiff=np.array(sorted([list(i.get()) for i in results],key=lambda x: x[0],reverse=False))
    
    opt_Imat = None
    
    if isinstance(optimal,int):
        ncomp_opt=optimal
    elif optimal==True:
        ncomp_opt=int(np.argmax(list_Idiff[:,1]))
        ncomp_opt=int(list_Idiff[ncomp_opt,0])
    
    if optimal != None:
        opt_Imat,_,_,_=compute_Imat_from_listFCs(fit_kcomponents(U,S,Vt,all_connectomes,mean_,ncomp_opt),comp_metric)
    
        
    return (list_Idiff), varexp, opt_Imat


def plot_classic_Imat_and_PCAcurve_concat(data1,data2,data,trimming_length=0,vmin=-0.2,vmax=1,condition='Rest',parallel=None,
                                          ticks_width=1.5, axis_width=2.5, labelsizereg=16,comp_metric="pearson_correlation"):
    '''
    Function for brain+spine appending the 2 matrixes discarding the between connections
    
    '''
    
    FC_list=computing_FCs_halved_concat(data1,data2,data,trimming_length=trimming_length)
    I_mat,Idiff,Iself,Iothers=compute_Imat_from_single_session_FCs(FC_list,comp_metric)
    cohend = cohen_d(I_mat.diagonal(),I_mat[~np.eye(I_mat.shape[0],dtype=bool)])
    
    fig = plt.figure(dpi=150,figsize=(12,6))
    ax= plt.subplot(121)
    
    im=plt.imshow(I_mat,vmin=vmin,vmax=vmax)
    cbar=plt.colorbar(im,fraction=0.046, pad=0.04)
    
    plt.xlabel('Subjects Test (day 1)')
    plt.ylabel('Subjects Retest (day 2)')
    cbar.set_label(r"Pearson's $\rho$", rotation=90,labelpad=-5)
    num_title=list([np.round(i,2) for i in [Idiff,cohend,Iself,Iothers]])
    plt.title('%s \n Idiff: %.2f %%, Cohen D: %.2f\nIself: %.2f, Iothers: %.2f' %(condition,num_title[0],num_title[1],num_title[2],num_title[3]));
    plot_adjustments(ax,ticks_width=ticks_width, axis_width=axis_width, labelsizereg=labelsizereg)
    ax.spines["top"].set_visible(True)
    ax.spines["right"].set_visible(True)


    ax= plt.subplot(122)

    PCA_list_idiff=compute_Idiff_PCA_one_dataset_from_FCs(FC_list)
    ax.plot(PCA_list_idiff[:,0],PCA_list_idiff[:,1],'ro-')
    asp = np.diff(ax.get_xlim())[0] / np.diff(ax.get_ylim())[0]
    ax.set_aspect(asp)
    plt.xlabel('Number of PCA components')
    plt.ylabel(r'$I_{diff}$ (%)',fontsize=14,labelpad=-1)
    plt.subplots_adjust(wspace=0.35)
    plot_adjustments(ax,ticks_width=ticks_width, axis_width=axis_width, labelsizereg=labelsizereg)
    return(I_mat,PCA_list_idiff)

    
    

def plot_classic_Imat_and_PCAcurve(data_test,data_retest=None,trimming_length=0,vmin=-0.2,vmax=1,condition='Rest',parallel=None,
                                   ticks_width=1.5, axis_width=2.5,labelsizereg=16,
                                   metric="pairwise_correlation",optimal=False, inter_bet_con=None,W=1,comp_metric="pearson_correlation"):
    '''
    Function that plot Imat and the Idiff values as a function of the number of components
    as proposed in the paper by Amico and Goni "The quest for identifiability
    in human functional connectomes".
    If only one dataset is present, then do the procedure when splitting the recordings in half

    input: data1 -> dictionary in the form {subjID: fMRI data} %% TEST
           data2 -> dictionary in the form {subjID: fMRI data} %% RETEST
           trimming_length -> default 0 (no trimming), if different from 0 th
    '''
    if data_retest==None:
        FC_list=computing_FCs_halved(data_test,trimming_length=trimming_length,metric=metric)   
        I_mat,Idiff,Iself,Iothers=compute_Imat_from_single_session_FCs(FC_list,comp_metric)
    else:
        FC1_list=computing_FCs(data_test,trimming_length=trimming_length,metric=metric,inter_bet_con=inter_bet_con, W=W)
        FC2_list=computing_FCs(data_retest,trimming_length=trimming_length,metric=metric,inter_bet_con=inter_bet_con, W=W)            
        I_mat,Idiff,Iself,Iothers=compute_Imat_from_FCs_twodatasets(FC1_list,FC2_list,comp_metric)

    cohend = cohen_d(I_mat.diagonal(),I_mat[~np.eye(I_mat.shape[0],dtype=bool)])
    fig = plt.figure(dpi=150,figsize=(12,6))
    ax= plt.subplot(121)
    
    im=plt.imshow(I_mat,vmin=vmin,vmax=vmax)
    cbar=plt.colorbar(im,fraction=0.046, pad=0.04)
    if data_retest==None:
        plt.xlabel('Subjects Test (first half)')
        plt.ylabel('Subjects Retest (second half)')
    else:
        plt.xlabel('Subjects Test (day 1)')
        plt.ylabel('Subjects Retest (day 2)')
    cbar.set_label(r"Pearson's $\rho$", rotation=90,labelpad=-5)
    num_title=list([np.round(i,2) for i in [Idiff,cohend,Iself,Iothers]])
    plt.title('%s \n Idiff: %.2f %%, Cohen D: %.2f \nIself: %.2f, Iothers: %.2f' %(condition,num_title[0],num_title[1],num_title[2],num_title[3]));
    plot_adjustments(ax,ticks_width=ticks_width, axis_width=axis_width, labelsizereg=labelsizereg)
    ax.spines["top"].set_visible(True)
    ax.spines["right"].set_visible(True)


    ax= plt.subplot(122)
    varexp = None
    if data_retest==None:
        PCA_list_idiff,Imat_opt=compute_Idiff_PCA_one_dataset_from_FCs(FC_list,optimal=optimal)
    else:
        PCA_list_idiff,varexp,Imat_opt=compute_Idiff_PCA_two_datasets_from_FCs(FC1_list,FC2_list,parallel=parallel,optimal=optimal)
        
    ax.plot(PCA_list_idiff[:,0],PCA_list_idiff[:,1],'ro-')
    asp = np.diff(ax.get_xlim())[0] / np.diff(ax.get_ylim())[0]
    ax.set_aspect(asp)
    plt.xlabel('Number of PCA components')
    plt.ylabel(r'$I_{diff}$ (%)',fontsize=14,labelpad=-1)
    plt.subplots_adjust(wspace=0.35)
    plot_adjustments(ax,ticks_width=ticks_width, axis_width=axis_width, labelsizereg=labelsizereg)
    
    return(I_mat,PCA_list_idiff,varexp, Imat_opt)




## ICC functions


def ICC(matrix, alpha=0.05, r0=0):
    '''Intraclass correlation
    matrix is matrix of observations. Each row is an object of measurement and
    each column is a judge or measurement.
    '1-1' is implement: The degree of absolute agreement among measurements made on
    randomly seleted objects. It estimates the correlation of any two
    measurements.
    ICC is the estimated intraclass correlation. LB and UB are upper
    and lower bounds of the ICC with alpha level of significance. 

    In addition to estimation of ICC, a hypothesis test is performed
    with the null hypothesis that ICC = r0. The F value, degrees of
    freedom and the corresponding p-value of the this test are reported.

    Translated in python from the matlab code of Arash Salarian, 2008

    Reference: McGraw, K. O., Wong, S. P., "Forming Inferences About
    Some Intraclass Correlation Coefficients", Psychological Methods,
    Vol. 1, No. 1, pp. 30-46, 1996'''

    M = np.array(matrix)
    n, k = np.shape(M)
    SStotal = np.var(M.flatten(), ddof=1) * (n * k - 1)
    MSR = np.var(np.mean(M, 1), ddof=1) * k
    MSW = np.sum(np.var(M, 1, ddof=1)) / n
    MSC = np.var(np.mean(M, 0), ddof=1) * n
    MSE = (SStotal - MSR * (n - 1) - MSC * (k - 1)) / ((n - 1) * (k - 1))

    r = (MSR - MSW) / (MSR + (k - 1) * MSW)

    F = (MSR / MSW) * (1 - r0) / (1 + (k - 1) * r0)
    df1 = n - 1
    df2 = n * (k - 1)

    p = 1 - stats.f.cdf(F, df1, df2)

    FL = (MSR / MSW) * (stats.f.isf(1 - alpha / 2, n * (k - 1), n - 1))
    FU = (MSR / MSW) / (stats.f.isf(1 - alpha / 2, n - 1, n * (k - 1)))

    LB = (FL - 1) / (FL + (k - 1))
    UB = (FU - 1) / (FU + (k - 1))

    return (r, LB, UB, F, df1, df2, p)
    

def compute_ICC_parallel(i,mat):
    return(i,ICC(mat)[0])
    

def compute_ICC_idenfiability_from_FCs(FC_lists,parallel=None):
    '''
    Compute the ICC matrix '1-1' considering the list of flattened functional connectomes
    '''
    
    num_FCs=len(FC_lists)
    nsubjs,edges=FC_lists[0].shape
    nROIs=int(np.ceil(np.sqrt(2*edges)))
    ICC_matrix=np.zeros((nROIs,nROIs))
    u,v= np.triu_indices(n=nROIs,k=1,m=nROIs)
    ICC_values=np.zeros((edges))
    if parallel==None:
        for i in range(edges):
            current_mat_edges=np.vstack([FC_lists[s][:,i] for s in range(num_FCs)])
            
            ICC_values[i]=ICC(current_mat_edges.T)[0]
    else:
        pool=Pool(processes=parallel)
        res=[]
        for i in range(edges):
            
            ICC_mat=np.vstack([FC_lists[s][:,i] for s in range(num_FCs)]).T
            res.append(
                pool.apply_async(compute_ICC_parallel,
                                 args=(i,ICC_mat)))
        pool.close()
        pool.join()
        list_ICC=np.array(sorted([list(i.get()) for i in res],key=lambda x: x[0],reverse=False))
#        print(list_ICC)
        ICC_values=list_ICC[:,1]
    ICC_matrix[u,v]=ICC_values
    ICC_matrix+=ICC_matrix.T
    return(ICC_matrix)
    
def compute_ICC_idenfiability_from_single_FCs(FC_list,parallel=None):
    '''
    Compute the ICC matrix '1-1' considering the list of flattened functional connectomes
    '''
    
    FC1_list=FC_list[:,0,:]
    FC2_list=FC_list[:,1,:]
    
    nsubjs,edges=FC1_list.shape
    nROIs=int(np.ceil(np.sqrt(2*edges)))
    ICC_matrix=np.zeros((nROIs,nROIs))
    u,v= np.triu_indices(n=nROIs,k=1,m=nROIs)
    
    ICC_values=np.zeros((edges))
    if parallel==None:
        for i in range(edges):
            ICC_values[i]=ICC(np.vstack((FC1_list[:,i],FC2_list[:,i])).T)[0]
    else:
        pool=Pool(processes=parallel)
        res=[]
        for i in range(edges):
            ICC_mat=np.vstack( (FC1_list[:,i],FC2_list[:,i]) ).T
            res.append(
                pool.apply_async(compute_ICC_parallel,
                                 args=(i,ICC_mat)))
        pool.close()
        pool.join()
        list_ICC=np.array(sorted([list(i.get()) for i in res],key=lambda x: x[0],reverse=False))
        ICC_values=list_ICC[:,1]
    ICC_matrix[u,v]=ICC_values
    ICC_matrix+=ICC_matrix.T
    return(ICC_matrix)




def find_limit_yeoOrder(yeovector,deflist=['VIS','SM','DA','VA','L','FP','DMN','C'],N=400,subcortical=None):
    scroll_deflist=0
    if subcortical!=None:
        list_yeo_ends={i:[0,0] for i in deflist}
    else:
        list_yeo_ends={i:[0,0] for i in deflist[:-1]}
    #print(list_yeo_ends)
    yeo_init=0
    for idx,i in enumerate(yeovector[:-1]):
        val_current=i
        val_prox=yeovector[idx+1]
        if val_prox-val_current < 0:
            yeo_end=idx
            list_yeo_ends[deflist[scroll_deflist]][0]=yeo_init
            list_yeo_ends[deflist[scroll_deflist]][1]=yeo_end+1
            scroll_deflist+=1
            yeo_init=idx+1
    ##This is for the DMN
    list_yeo_ends[deflist[scroll_deflist]][0]=yeo_init
    list_yeo_ends[deflist[scroll_deflist]][1]=N
    if subcortical != None:
        list_yeo_ends[deflist[scroll_deflist+1]][0]=N+1
        list_yeo_ends[deflist[scroll_deflist+1]][1]=len(yeovector)
    return(list_yeo_ends) 



def load_yeonets_matfile(filename='/home/andrea/Dropbox/Brain_data/Misc/yeoOrder/yeo_RS7_Schaefer400S.mat',subcortical=None):
    Yeonets=loadmat(filename)
    N_yeos=len(Yeonets['yeoOrder'])
    if subcortical==None:
        yeoOrder=np.array([i[0]-1 for i in Yeonets['yeoOrder']])[:int(N_yeos-19)]
        limit_yeo=find_limit_yeoOrder(yeoOrder,N=int(N_yeos-19),subcortical=None)
    else:
        yeoOrder=np.array([i[0]-1 for i in Yeonets['yeoOrder']])
        limit_yeo=find_limit_yeoOrder(yeoOrder,N=N_yeos,subcortical=True)
    return(yeoOrder,limit_yeo)


def plot_ICC_mat(ICC_mat,yeo_net=True,yeoOrder=None,limit_yeo=None,
                subcortical=None,percentile=95,ax=None):
    if ax == None:
        fig=plt.figure(dpi=150)
        ax=plt.subplot(111)
    
    list_colors_yeo_center=['#d9d9d9','#7fc97f', '#beaed4',
                     '#fdc086', '#ffeda0', '#386cb0',
                     '#f0027f','#bf5b17','#d9d910']
    list_colors_yeo=['#969696','#7fc97f', '#beaed4',
                     '#fdc086','#ffeda0', '#386cb0',
                     '#f0027f','#bf5b17']
    cmap = colors.ListedColormap(list_colors_yeo)

    if yeo_net==None:
        cmap =  colors.ListedColormap(['darkgray','white'])
        cleaned_ICC_mat=np.where(ICC_mat>np.percentile(ICC_mat,q=percentile),1,np.nan)
        ax.matshow(cleaned_ICC_mat,cmap=cmap)
    else:
        ICC_mat_reodered=ICC_mat[yeoOrder,:][:,yeoOrder]
        yeo_ICC_mat=np.where(ICC_mat_reodered>np.percentile(ICC_mat,q=percentile),0,np.nan)
        k=1
        for net,limits in limit_yeo.items():
            y0,y1=limits[0],limits[1]
#             print(y0,y1,k)
            yeo_ICC_mat[y0:y1,y0:y1]=np.where(yeo_ICC_mat[y0:y1,y0:y1]==0,k,np.nan)
            k+=1

        ax.matshow(yeo_ICC_mat,cmap=cmap)
        ax.set_xlim(-0,len(ICC_mat));
        ax.set_ylim(len(ICC_mat),0);
        s=0
        net_names=['VIS','SM','DA','VA','L','FP','DMN','SC']
        for net,limits in limit_yeo.items():
            y0,y1=limits[0],limits[1]
#             if y0 <= len(ICC_mat) and y1<= len(ICC_mat):
                #print(y0,y1)
            rect_upper = Rectangle((len(ICC_mat)+2,y0),5, y1-y0, color =list_colors_yeo_center[s+1], clip_on=False,zorder=10)
            rect_bottom = Rectangle((y0,len(ICC_mat)+2),y1-y0,5,color =list_colors_yeo_center[s+1], clip_on=False,zorder=10)
            ax.text(x=1.025*len(ICC_mat),y=y0+(0.65*(y1-y0)),s=net_names[s],clip_on=False)
            s+=1
            ax.add_patch(rect_upper)
            ax.add_patch(rect_bottom)
            plt.hlines(y0,y0,y1,colors='k',lw=1)
            plt.hlines(y1,y0,y1,colors='k',lw=1)
            plt.vlines(y0,y0,y1,colors='k',lw=1)
            plt.vlines(y1,y0,y1,colors='k',lw=1)



def compute_ICC_yeo_net_dataframe(ICC_mat,yeoOrder,limit_yeo):
    ICC_mat_reodered=ICC_mat[yeoOrder,:][:,yeoOrder]
    list_average_ICC=[]
    list_yeo_repeated=[]
    s=1
    for net,limits in limit_yeo.items():
        y0,y1=limits[0],limits[1]
        gap=y1-y0
#         if y0 <= len(ICC_mat) and y1<= len(ICC_mat):
#         print(y0,y1,len(ICC_mat_reodered[y0:y1,y0:y1][~np.eye(gap,dtype=bool)].flatten()),(gap)*(gap-1))
        list_average_ICC.extend(ICC_mat_reodered[y0:y1,y0:y1][~np.eye(gap,dtype=bool)].flatten())
        list_yeo_repeated.extend(np.repeat(s,(gap)*(gap-1)))
        s+=1
    current_df=pd.DataFrame([list_average_ICC,list_yeo_repeated]).T
    current_df=current_df.rename(columns={0:'ICC',1:'YeoNet'})
    return(current_df)

def plot_ICC_violins(ICC_mat,yeoOrder,limit_yeo,ax=None):
    current_df=compute_ICC_yeo_net_dataframe(ICC_mat,yeoOrder,limit_yeo)
    list_colors_yeo=['#969696','#7fc97f', '#beaed4',
                     '#fdc086','#ffeda0', '#386cb0',
                     '#f0027f','#bf5b17']
    if ax == None:
        fig=plt.figure(dpi=150)
        ax=plt.subplot(111)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

    sns.violinplot(x=current_df['YeoNet'],y=current_df['ICC'],orient='v',cut=0.5,
                   palette=sns.color_palette(list_colors_yeo[1:]),alpha=0.5,width=0.7,linewidth=1)
    
    net_names=list(limit_yeo.keys())
    # Set the x labels
    ax.set_ylabel('ICC Values',fontsize=14)
    ax.set_xlabel('')
    ax.set_ylim(0,1)
    ax.set_xticklabels(net_names);


def plot_ICC_mat_and_violins(ICC_mat,yeoOrder,limit_yeo,percentile=95,subcortical=None):
    fig=plt.figure(figsize=(10,6),dpi=150)
    gs = GridSpec(6, 2)
    ax=fig.add_subplot(gs[0:5, 0])
    plot_ICC_mat(ICC_mat,yeo_net=True,yeoOrder=yeoOrder,subcortical=subcortical,
                limit_yeo=limit_yeo,percentile=percentile,ax=ax)

    ax=fig.add_subplot(gs[1:4, 1])
    plot_ICC_violins(ICC_mat,yeoOrder=yeoOrder,limit_yeo=limit_yeo,ax=ax)
    plt.subplots_adjust(wspace=0.32)