From eff7e8bc577a9e86e8ad1bbcb5c46a06a7a0338f Mon Sep 17 00:00:00 2001 From: Claire Shao Date: Wed, 28 May 2025 05:11:50 -0400 Subject: [PATCH] Add deep learning 2022 examples --- examples/eg__griffiths2022.py | 2064 +++++++++++++++++++++++++++++++++ 1 file changed, 2064 insertions(+) create mode 100644 examples/eg__griffiths2022.py diff --git a/examples/eg__griffiths2022.py b/examples/eg__griffiths2022.py new file mode 100644 index 00000000..f06e1c52 --- /dev/null +++ b/examples/eg__griffiths2022.py @@ -0,0 +1,2064 @@ +# -*- coding: utf-8 -*- +"""eg__griffiths2022.ipynb + +Automatically generated by Colab. + +Original file is located at + https://colab.research.google.com/drive/1ijRJcgcrtjDpP1r9XyUCVLT19E8hUPwC + +#0. Overview +This example contains the replication results of the paper: + Griffiths, J. D., Wang, Z., Ather, S. H., Momi, D., Rich, S., Diaconescu, A., McIntosh, A. R., & Shen, K. + Deep Learning-Based Parameter Estimation for Neurophysiological Models of Neuroimaging Data. + bioRxiv (2022). + DOI: 10.1101/2022.05.19.492664 + + +The code in this example allows for the replication of the findings in the paper. If you utilize this code, please cite the original publication. + +# 1. Setup +General and whobpyt importage +""" + +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from scipy.optimize import fsolve +import sys +import json +import pandas as pd # for data manipulation +import seaborn as sns # for plotting +import time # for timer +import os +from numpy import exp,sin,cos,sqrt,pi, r_,floor,zeros +import scipy.io +from matplotlib.tri import Triangulation +#from nilearn import plotting, datasets +from matplotlib import cm +from matplotlib.pyplot import subplot +import warnings # for suppressing warnings and output +import nibabel as nib +import matplotlib.pyplot as plt # for plotting +warnings.filterwarnings('ignore') + +#sys.path.append('/content/drive/MyDrive/ClaireShao_WhoBPyT_Replications_Project/Paper 3 - dl-paramest-for-neurophys-models/dl-2024') + +from ≈.depr.griffiths2022.plotting import ( + get_rotation_matrix, + get_combined_rotation_matrix, + plot_surface_mpl_mv, + plot_surface_mpl, + plot_sim_states_outputs, + plot_fit_parameters, + R_stat +) +from whobpyt.depr.griffiths2022.bifurcation_analysis import ( + h_tf, + dh_tf, + smooth_normalize, + derivative_orig, + derivative, + get_eig_sys, + regime_search_I0, + regime_search_gEE +) +from whobpyt.depr.griffiths2022.rww_pytorch_model_updated import ( + RNNWWD, + Model_fitting +) +from whobpyt.depr.griffiths2022.wwd_test import WWD_test + +"""# 2. bifurcation_analysis_singlenode""" + +param = { + +# Parameters for the integration +"dt" : 0.001, # integration step size +"ROI_num" : 1, # number of neural nodes + + +# Parameters for the ODEs +# Excitatory population +"W_E" : 1., # scale of the external input +"tau_E" : 100., # decay time +"gamma_E" : 0.641/1000., # other dynamic parameter (?) + +# Inhibitory population +"W_I" : 0.7, # scale of the external input +"tau_I" : 10., # decay time +"gamma_I" : 1./1000., # other dynamic parameter (?) + +# External input +"I_0" : 0.45, # external input +"I_external" : 0., # external stimulation + +# Coupling parameters +"g" : 20., # global coupling (from all nodes E_j to single node E_i) +"g_EE" : .5, # local self excitatory feedback (from E_i to E_i) +"g_IE" : .14, # local inhibitory coupling (from I_i to E_i) +"g_EI" : 0.15, # local excitatory coupling (from E_i to I_i) +"lamb" : 0., # scale of global coupling on I_i (compared to E_i) + +"aE":310, +"bE" :125, +"dE":0.16, +"aI":615, +"bI" :177, +"dI" :0.087, + + + +} ### initial values of all model prameters + +# strategies +# Brutal search: +#firstly" find largest I_0 that the system has not bistable regym +### +"""param_study = ['g_EE', 'g_IE', 'g_EI', 'I_0'] +num_param = 4 +num_trials = 20 + +I0_rng = np.linspace(0.35,.1, 5) +gEE_rng = np.linspace(0.4,1.5, num_trials) +gIE_rng = np.linspace(0.1,0.5, num_trials) +gEI_rng = np.linspace(0.1,0.5, num_trials) +c_I0 = regime_search_I0(I0_rng, gEE_rng, gIE_rng, gEI_rng, param) #CS +I0 = I0_rng[np.argmax(c_I0)] +print('I0 <', I0, 'bifurcation') +""" + +file_path ='/content/drive/MyDrive/ClaireShao_WhoBPyT_Replications_Project/Paper 3 - dl-paramest-for-neurophys-models/regime_search_I0_result.json' +with open(file_path, 'r') as f: + data = json.load(f) + +# Now access I0 and c_I0 +I0 = data['I0'] +c_I0 = data['c_I0'] +print('I0 <', I0, 'bifurcation') + +# fixed I_0 = 0.2 coastly search g_EE + +param_study = ['g_EE', 'g_IE', 'g_EI', 'I_0'] +num_param = 4 +num_trials = 20 +I0_rng = np.linspace(0.2,.25, 1) +gEE_rng = np.linspace(0.4,1, 10) +gIE_rng = np.linspace(0.1,0.5, num_trials) +gEI_rng = np.linspace(0.1,0.5, num_trials) + +c_gEE = regime_search_gEE(I0_rng, gEE_rng, gIE_rng, gEI_rng, param) #CS +gEE = gEE_rng[np.argmax(c_gEE)] +print('gEE >', gEE, 'bistable') + +#### I_0 = 0.2 gEE=0.5 +# seach gEI and gIE + +param_study = ['g_EE', 'g_IE', 'g_EI']#, 'I_0'] +num_param = 2 +num_trials = 50 + + +#gEE_rng = np.linspace(0.2,1, 20) +gIE_rng = np.linspace(0.0,2, num_trials) +gEI_rng = np.linspace(0.0,2, num_trials) +fig = plt.figure() +ax = fig.add_subplot(111)#, projection='3d') + +gIE, gEI = np.meshgrid(gIE_rng, gEI_rng)#, I0_rng) +#gIE, gEI = np.meshgrid(gIE_rng, gEI_rng)#, I0_rng) +gEE_good =[] +gIE_good =[] +gEI_good =[] +#I0_good =[] +c =[] +eig_v=[] + + +param['I_0'] = 0.2 +param['g_EE'] = 0.5 + + +for i in range(num_trials**num_param): + ind_0 = i//num_trials**(num_param -1) + ind_1 = (i%num_trials**(num_param -1)) //num_trials**(num_param -2) + #ind_2 = ((i%num_trials**(num_param -1)) % num_trials**(num_param -2)) // num_trials**(num_param -3) + #ind_3 = ((i%num_trials**(num_param -1)) % num_trials**(num_param -2)) % num_trials**(num_param -3) + #param['g_EE'] = gEE + param['g_IE'] = gIE[ind_0][ind_1]#[ind_2]#[ind_3] + param['g_EI'] = gEI[ind_0][ind_1]#[ind_2]#[ind_3] + #param['I_0'] = I0[ind_0][ind_1][ind_2][ind_3] + #param['g_EE'] = gEE[ind_0][ind_1][ind_2]#[ind_3] + """param['g_IE'] = gIE[ind_0][ind_1]#[ind_2]#[ind_3] + param['g_EI'] = gEI[ind_0][ind_1]#[ind_2]#[ind_3]""" + + """gEE_good.append(param['g_EE'])""" + gIE_good.append(param['g_IE']) + gEI_good.append(param['g_EI']) + #print(param['g_EE'], param['g_IE'], param['g_EI'], param['I_0']) + initial = np.random.uniform(0.,2, [2, 200]) + solns = [] + for i in range(initial.shape[1]): + x0 = initial[:,i] + #print('x0', x0) + for j in range(20): + x0 = np.round(fsolve(derivative, x0, param), decimals =4) + + #print('solutions', solutions) + #print(derivative(solutions)) + if (np.abs(derivative(x0, param)) > 1.0).sum() == 0: + solns.append(tuple(x0)) + #print(len(set(solns))) + good_sols =[] + for sol in set(solns): + sol_good = True + if len(good_sols)>0: + for g_sol in good_sols: + if np.sqrt(((np.array(g_sol)-np.array(sol))**2).mean()) < 1e-3: + sol_good = False + break + if sol_good == True: + good_sols.append(sol) + + #print(len(good_sols)) + c.append(len(good_sols)) + eig_ins=[] + for sol in good_sols: + #print(sol, derivative(np.array(sol))) + eig = get_eig_sys(sol[0], sol[1], param) + eig_ins.append(eig) + eig_v.append(eig_ins) +im = ax.scatter(np.array(gIE_good), np.array(gEI_good), cmap ='RdYlBu_r', c = np.array(c), vmax=4, vmin=0) +fig.colorbar(im, ax = ax) +plt.show() + +data = np.load('bifurcation_results.npz', allow_pickle=True) + +gIE_good = data['gIE_good'] +gEI_good = data['gEI_good'] +c = data['c'] +eig_v = data['eig_v'] # only if needed + +# Plot again +fig = plt.figure() +ax = fig.add_subplot(111) +im = ax.scatter(gIE_good, gEI_good, cmap='RdYlBu_r', c=c, vmax=4, vmin=0) +fig.colorbar(im, ax=ax) +plt.title('Bifurcation Diagram (from saved result)') +plt.show() + +# Figures for the paper +fig = plt.figure(figsize=(6,4)) +ax = fig.add_subplot(111) +im = ax.scatter(np.array(gIE_good), np.array(gEI_good), cmap ='bwr', c = np.array(c), vmax =2, vmin=1, alpha=0.3) +ax.set_xlabel('$g_{IE}$', fontsize=12) +ax.set_ylabel('$g_{EI}$', fontsize=12) +ax.text(1,1.25, 'I', fontsize=15, weight='bold') +ax.text(.1,.15, 'II', fontsize=15, weight='bold') +#fig.colorbar(im, ax = ax) +ax.set_title("Bifurcation analysis on $g_{IE}$ and $g_{EI}$ when $g_{EE} = 0.5$ and $I_0=0.2$") +#fig.savefig('/home/jwang/bifurcationanalysis_gEI_gIE.png') +plt.show() + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +I0= 0.4167 +x1=np.linspace(0,I0- 0.001,50) +x2 =np.linspace(I0,2) +legend_text =['stable fixed points', 'unstable fixed points'] +y1 = np.sqrt(-x1+I0) +y2= - np.sqrt(-x1+I0) +ax.plot(x1, y1, 'r', linewidth=10) +ax.plot(x1, y2, 'r', linewidth=10) +h1, = ax.plot(x1, np.zeros_like(x1), 'r', linewidth=10) +h2, = ax.plot(x2,np.zeros_like(x2), 'g', linewidth=10) +ax.text(I0+0.05, 0.15, '$I_0=$'+str(I0), fontsize = 30) +ax.tick_params(labelsize = 30) +ax.legend([h2,h1], legend_text, fontsize = 20) +plt.show() + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +gEE= 0.4526 +x1=np.linspace(0,gEE,50) +x2 =np.linspace(gEE+0.001,2) +legend_text =['stable fixed points', 'unstable fixed points'] +y1 = np.sqrt(x2-gEE) +y2= - np.sqrt(x2-gEE) +ax.plot(x2, y1, 'g', linewidth=10) +ax.plot(x2, y2, 'g', linewidth=10) +h1, = ax.plot(x1, np.zeros_like(x1), 'g', linewidth=10) +h2, = ax.plot(x2,np.zeros_like(x2), 'r', linewidth=10) +ax.text(gEE+0.05, 0.15, '$g_{EE}=$'+str(gEE), fontsize = 30) +ax.tick_params(labelsize = 30) +ax.legend([h1,h2], legend_text, fontsize = 20) +plt.show() + +"""# 3. Synthetic data model-fitting analysis +The initial sections of the code handle downloading and importing the required files. Please ensure you manually download the necessary files and update the file paths accordingly. + +Additionally, for compatibility with newer versions of NumPy, if you encounter errors related to numpy.int, please update it to numpy.int64 in rww_pytorch_model.py + +### 3a +""" + +mask = np.tril_indices(66,-1) + +node_size = 66 +output_size = 66 +num_epoches = 1 +batch_size = 20 +step_size = .05 +input_size = 2 +tr = .75 + +G = 100 +gEE = 3.5 +gIE = 0.42 +gEI = 0.42 + + +sc_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/weights.txt' #CS +sc = np.loadtxt(sc_file) +sc = sc -np.diag(np.diag(sc)) +sc = 0.5*(sc.T+sc) +sc = np.log1p(sc)/np.linalg.norm(np.log1p(sc)) + +# make a fake ts +len_sim = 2400 +ts = np.ones((len_sim, node_size)) + + + +model = RNNWWD(input_size, node_size, batch_size, step_size, tr, sc, False, g_mean_ini=100, g_std_ini = .1, gEE_mean_ini=2.5, gEE_std_ini = .1) + + +# call model fit method +F = Model_fitting(model, ts, num_epoches) + + +output_test = F.test(240) + +fc_sim = np.corrcoef(output_test['simBOLD'][:,2400:]) + +## connectivity 66 +fig, ax=plt.subplots(1,2, figsize=(12,8)) +ax[0].imshow(fc_sim - np.diag(np.diag(fc_sim)), cmap='bwr') +ax[1].imshow(sc, cmap='bwr') +plt.show() + +"""###3b""" + +out_dir = 'results/' +if not os.path.isdir(out_dir): + os.makedirs(out_dir) + +base_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP' #CS + +subs =sorted([sc_file[-10:-4] for sc_file in os.listdir(base_dir) if sc_file[:8] == 'weights_']) + +start_time = time.time() + +#subs_s = ['562446', '257542', '154936' ] +base_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/' +out_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/' +for i in range(0,40): + + + node_size = 83 + mask = np.tril_indices(node_size, -1) + num_epoches = 20 + batch_size = 20 + step_size = 0.05 + input_size = 2 + tr = 0.75 + sub=subs[i] + print(i, sub) + sc_file = base_dir +'weights_'+sub+'.txt' + ts_file = base_dir +sub+'_rfMRI_REST1_LR_hpc200_clean__l2k8_sc33_ts.pkl'#out_dir+'sub_'+sub+'simBOLD_idt.txt'# + + if os.path.isfile(sc_file) and os.path.isfile(ts_file): + sc = np.loadtxt(sc_file) + SC =(sc+sc.T)*0.5 + + sc = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) + + ts_pd = pd.read_pickle(ts_file) + ts = ts_pd.values + #ts = np.loadtxt(ts_file) + ts =ts/np.max(ts) + fc_emp = np.corrcoef(ts.T) + # Get the WWD model module for forward in a batch. + + + model = RNNWWD(input_size, node_size, batch_size, step_size, tr, sc, True, g_mean_ini=80, g_std_ini = .1, gEE_mean_ini=2.5, gEE_std_ini = .1) + + + + # call model fit method + F = Model_fitting(model, ts, num_epoches) + + # fit data(train) + + output_train = F.train() + + + output_test = F.test(120) + plot_fit_parameters(output_train) + plot_sim_states_outputs(ts, output_test) + + sc_mod = np.zeros_like(sc) + sc_mod[mask] = output_train['gains'][-100:].mean(0) + sc_mod = sc_mod+sc_mod.T + w = (1 + np.tanh(sc_mod))*sc + w_n = 0.5*(w + w.T)/np.linalg.norm(0.5*(w + w.T)) + np.savetxt(out_dir + 'bold_test_'+ sub +'.txt', output_test['simBOLD']) + np.savetxt(out_dir + 'bold_train_'+ sub +'.txt', output_train['simBOLD']) + np.savetxt(out_dir + 'sc_mod_'+ sub +'.txt', sc_mod) + np.savetxt(out_dir + 'sc_'+ sub +'.txt', sc) + g= output_train['g'][-100:].mean() + gEE = output_train['gEE'][-100:].mean() + gEI = output_train['gEI'][-100:].mean() + gIE = output_train['gIE'][-100:].mean() + + + np.savetxt(out_dir + 'parameters_'+ sub +'.txt', np.array([g,gEE, gIE, gEI]).T) +end_time = time.time() +print('running time is {0} \'s'.format(end_time - start_time )) + +Syn_ts_sim ={} +Syn_ts_test ={} +Syn_paras = {} +Syn_sc_mod = {} +Syn_sc = {} +for i in range(39): + sub = 'sub_'+str(i) + + + Syn_ts_test[sub] = np.loadtxt('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/bold_test_'+ sub +'.txt') + Syn_ts_sim[sub] = np.loadtxt('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann//bold_train_'+ sub +'.txt') + Syn_sc_mod[sub] = np.loadtxt('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/sc_mod_'+ sub +'.txt') + Syn_sc[sub] = np.loadtxt('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann//Hagmann/sc_'+ sub +'.txt') + + Syn_paras[sub] = np.loadtxt('../data/Hagmann/parameters_'+ sub +'.txt') + +np.save('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/Syn_ts_sim.npy', Syn_ts_sim) +np.save('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/Syn_ts_test.npy', Syn_ts_test) +np.save('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/Syn_sc_mod.npy', Syn_sc_mod) +np.save('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/Syn_sc.npy', Syn_sc) +np.save('/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/Syn_fitparas.npy', Syn_paras) + +"""### 3c""" + +mask = np.tril_indices(66,-1) + +node_size = 83 +step_size = 0.05 +Tr = 2.5 + +g = 82 +gEE = 2.5 +gIE = 0.4473 +gEI = 0.501 + +#CS here changed the direction of her own, it might be different than others +sc_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/weights.txt' +sc = np.loadtxt(sc_file) +sc = (sc + sc.T)*0.5 +SC = sc - np.diag(np.diag(sc)) +#SC[SC 0].mean() +1.1*SC[SC>0].std()] =0 # for connectivity 96 only +wo = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) + +test_sim = WWD_test(g, gEE, gIE, gEI, wo, step_size, node_size, Tr) + +bold_test = test_sim.generate_bold(4800) #could be any number you want + +fc_sim = np.corrcoef(bold_test.T) + +## connectivity 83 +fig, ax=plt.subplots(1,2, figsize=(12,8)) +ax[0].imshow(fc_sim - np.diag(np.diag(fc_sim)), cmap='bwr') +ax[1].imshow(wo, cmap='bwr') +plt.show() + +ts_file = "SYN_bold_test_output.txt" # Change the name if needed + +# Save the array to a text file +np.savetxt(ts_file, bold_test) + +print(f"Saved bold_test results to {ts_file}") + +base_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/' #cs +#ts_file = base_dir+'syn_ts.txt' CS +ts_file = "/content/SYN_bold_test_output.txt" +#SC = np.loadtxt(SC_file) +ts = np.loadtxt(ts_file) +print(ts.shape) +sc_file = base_dir + 'weights.txt' + +sc = np.loadtxt(sc_file) +sc = (sc + sc.T)*0.5 +SC = sc - np.diag(np.diag(sc)) +#SC[SC 0].mean() + 1.*SC[SC>0].std()] = 0 +sc_true = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) #SC/np.linalg.norm(SC) +fc_true = np.corrcoef(ts[2400:3600,:].T) +fc_cross = np.corrcoef(ts[0:2400,:].T) + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +im = ax.imshow(sc_true, cmap='bwr') +mn = round(sc_true.min(),2) +mx = round(sc_true.max(), 2) +md = round((mx + mn)/2, 2) +mnd = round((md + mn)/2, 2) +mxd = round((mx + md)/2, 2) +cbar = plt.colorbar(im, ax=ax) +cbar.ax.tick_params(labelsize=30) +cbar.set_ticks([mn,mnd, md,mxd, mx]) +cbar.set_ticklabels([mn,mnd, md,mxd, mx]) +ax.set_xticklabels([]) +ax.set_yticklabels([]) +plt.show() +# plt.savefig(base_dir+'Syn_sc_true.png') + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +fc_nd = fc_true - np.diag(np.diag(fc_true)) +im = ax.imshow(fc_nd, cmap='bwr') +mn = round(fc_nd.min(),2) +mx = round(fc_nd.max(), 2) +md = round((mx + mn)/2, 2) +mnd = round((md + mn)/2, 2) +mxd = round((mx + md)/2, 2) +cbar = plt.colorbar(im, ax=ax) +cbar.ax.tick_params(labelsize=30) +cbar.set_ticks([mn,mnd, md,mxd, mx]) +cbar.set_ticklabels([mn,mnd, md,mxd, mx]) +ax.set_xticklabels([]) +ax.set_yticklabels([]) +plt.show() +# plt.savefig(base_dir+'Syn_fc_true.png') + +### Read the model-fitting data with gains +base_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/Hagmann/' +ts_sim_file = base_dir + 'Syn_ts_sim.npy' +ts_test_file = base_dir + 'Syn_ts_test.npy' +sc_mod_file = base_dir + 'Syn_sc_mod.npy' +sc_file = base_dir + 'Syn_sc.npy' +paras_file = base_dir + 'Syn_fitparas.npy' + +Syn_ts_sim_data = np.load(ts_sim_file, allow_pickle=True) +Syn_ts_test_data = np.load(ts_test_file, allow_pickle=True) +Syn_sc_mod_data = np.load(sc_mod_file, allow_pickle=True) +Syn_sc_data = np.load(sc_file, allow_pickle=True) +Syn_para_data = np.load(paras_file, allow_pickle=True) + +Syn_ts_test = Syn_ts_test_data.item() +Syn_ts_sim = Syn_ts_sim_data.item() +Syn_sc_mod = Syn_sc_mod_data.item() +Syn_sc = Syn_sc_data.item() +Syn_paras = Syn_para_data.item() + +start_time = time.time() + +plt.rcParams["axes.labelsize"] = 30 +mask = np.tril_indices(66,-1) +data_dict = {} +data_max = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] + +data_max['g'] = [] +data_max['gEE'] = [] +data_max['gIE'] = [] +data_max['gEI'] = [] + +""" +Across the test subjects, take the maximum of the absolute value of the mean of the parameter +values identified for the model and the same for the parameter values. Append the two values and +add them to the `data_max` dictionary for the coupling strength parameters. +""" + +para_name = ['g', + 'gEE', + 'gIE', + 'gEI'] +for sub in Syn_paras: + + data_max['g'].append(max(np.abs(Syn_paras[sub][0].mean()), 100)) + data_max['gEE'].append(max(np.abs(Syn_paras[sub][1].mean()), 3.5)) + data_max['gIE'].append(max(np.abs(Syn_paras[sub][2].mean()), 0.42)) + data_max['gEI'].append(max(np.abs(Syn_paras[sub][3].mean()), 0.42)) + +""" +Take these parameter values and use these maximum values to normalize the range of parameter values. +""" + +for sub in Syn_paras: + + #data_dict['g'].append((Syn_paras[sub][0].mean() - 100)/max(data_max['g'])) + data_dict['gEE'].append((Syn_paras[sub][1].mean() - 23)/max(data_max['gEE'])) #CS + data_dict['gIE'].append((Syn_paras[sub][2].mean() - 18)/max(data_max['gIE']))#CS + data_dict['gEI'].append((Syn_paras[sub][3].mean() - 17)/max(data_max['gEI']))#CS + data_dict['g'].append((Syn_paras[sub][0].mean() - 20)/max(data_max['g'])) + +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) + +""" +Plot the results. +""" +import seaborn as sns +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +# Ensure ax is properly indexed +fig, ax = plt.subplots(2, 2, figsize=(12, 8), sharex=True) +ax = np.array(ax) # Ensure proper indexing +plt.rcParams["axes.labelsize"] = 30 + +# Verify data range and adjust if needed +for i, key in enumerate(para_name): + data = np.array(data_dict[key]) + + if len(data) > 1: # KDE requires at least two points + print(f"Plotting {key}: min={data.min()}, max={data.max()}, std={data.std()}") # Debug info + + # Use a larger bandwidth for smoother KDE + sns.kdeplot(ax=ax[i // 2, i % 2], data=data, color='r', linewidth=3, bw_adjust=0.5) + + # Fallback: if KDE fails, show histogram + if not np.any(data): # Check if data is effectively empty + print(f"Warning: {key} data has no variance, switching to histogram.") + ax[i // 2, i % 2].hist(data, bins=20, color='r', alpha=0.7) + + ax[i // 2, i % 2].set_xlim(-0.3, 0.5) # Adjust x limits for visibility + ax[i // 2, i % 2].tick_params(labelsize=20) + ax[i // 2, i % 2].grid() + else: + print(f"Skipping {key} due to insufficient data for KDE.") + +fig.tight_layout() +plt.show() + +# Second Plot +fig, ax = plt.subplots(figsize=(12, 8)) +plt.rcParams["axes.labelsize"] = 25 + +if any(len(data_dict[key]) > 1 for key in data_dict): + sns.kdeplot(data=pd.DataFrame(data_dict), linewidth=4, bw_adjust=0.5) + ax.tick_params(labelsize=20) + ax.set_xlim(-1, 1) # Adjust x limits for visibility + ax.set_ylim(0, 20) # Adjust x limits for visibility + ax.grid() + plt.setp(ax.get_legend().get_texts(), fontsize='20') # Reduce fontsize for better visibility + plt.setp(ax.get_legend().get_title(), fontsize='25') + plt.show() + +start_time = time.time() +# plot hist of r: model fitting with gains +mask = np.tril_indices(66,-1) +corr_sim ={} +r =[] +r1 =[] +for sub in Syn_ts_sim: + r.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_sim[sub])[mask])[0,1]) + r1.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_test[sub])[mask])[0,1]) +corr_sim['train'] = np.array(r) +corr_sim['test'] = np.array(r1) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 35 +sns.kdeplot(data = corr_sim, linewidth = 10) +#ax.set_title('hist of Pearlson correlation between the fitting and empirical BOLD', fontsize =20) +ax.tick_params(labelsize= 25) +ax.grid() + +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize ='50') + +start_time = time.time() +# plot hist of r: model fitting with gains +mask = np.tril_indices(66,-1) +corr_sim ={} +r =[] +r1 =[] +for sub in Syn_ts_sim: + r.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_sim[sub][:,10:])[mask])[0,1]) + r1.append(np.corrcoef(fc_cross[mask], np.corrcoef(Syn_ts_test[sub])[mask])[0,1]) +corr_sim['train'] = np.array(r) +corr_sim['test'] = np.array(r1) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 35 +sns.kdeplot(data = corr_sim, linewidth = 10) +#ax.set_title('hist of Pearlson correlation between the fitting and empirical BOLD', fontsize =20) +ax.tick_params(labelsize= 25) +ax.grid() + +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize ='50') + +start_time = time.time() +# plot hist of r: model fitting with gains +mask = np.tril_indices(66,-1) +corr_sim ={} +r =[] +r1 =[] +for sub in Syn_ts_sim: + r.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_sim[sub][:,10:])[mask])[0,1]) + r1.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_test[sub])[mask])[0,1]) +corr_sim['train'] = np.array(r) +corr_sim['test'] = np.array(r1) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 35 +sns.kdeplot(data = corr_sim, linewidth = 10) +#ax.set_title('hist of Pearlson correlation between the fitting and empirical BOLD', fontsize =20) +ax.tick_params(labelsize= 25) +ax.grid() + +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize ='50') + +# plt.savefig(base_dir+'Syn_hist_fittingcorrelation.png') + +start_time = time.time() +mask = np.tril_indices(66,-1) +subs = ['randSClmtd', + 'diagrandSClmtd', + 'randSC', + 'dist'] + +titles = ['RandNonzeros', + 'RandNonzeroes+diag', + 'Rand', + 'DistRecp'] + +sc_diff_file = base_dir + 'Syn_diff_sc.npy' +fc_diff_sim_file = base_dir + 'Syn_diff_fc_sim.npy' +sc_diff_mod_file = base_dir + 'Syn_diff_sc_mod.npy' + +Syn_diff_sc_data = np.load(sc_diff_file, allow_pickle=True) +Syn_diff_fc_sim_data = np.load(fc_diff_sim_file, allow_pickle=True) +Syn_diff_sc_mod_data = np.load(sc_diff_mod_file, allow_pickle=True) +Syn_diff_sc = Syn_diff_sc_data.item() +Syn_diff_fc_sim = Syn_diff_fc_sim_data.item() + +Syn_diff_sc_mod = Syn_diff_sc_mod_data.item() +vmax_ls = [0,0,0] +vmin_ls = [1,1,1] +delta = np.array([0,0,0,0,-0.001]) +i = 0 +for sub in Syn_diff_sc_mod: + if i < 4: + + if Syn_diff_sc[sub].max() > vmax_ls[0]: + vmax_ls[0] = Syn_diff_sc[sub].max() + if Syn_diff_sc_mod[sub].max() > vmax_ls[1]: + vmax_ls[1] = Syn_diff_sc_mod[sub].max() + if (Syn_diff_fc_sim[sub]- np.diag(np.diag(Syn_diff_fc_sim[sub]))).max() > vmax_ls[2]: + vmax_ls[2] = (Syn_diff_fc_sim[sub]- np.diag(np.diag(Syn_diff_fc_sim[sub]))).max() + + if Syn_diff_sc[sub].min() < vmin_ls[0]: + vmin_ls[0] = Syn_diff_sc[sub].min() + if Syn_diff_sc_mod[sub].min() < vmin_ls[1]: + vmin_ls[1] = Syn_diff_sc_mod[sub].min() + if (Syn_diff_fc_sim[sub]- np.diag(np.diag(Syn_diff_fc_sim[sub]))).min() < vmin_ls[2]: + vmin_ls[2] = (Syn_diff_fc_sim[sub]- np.diag(np.diag(Syn_diff_fc_sim[sub]))).min() + i += 1 +fig, ax = plt.subplots(4, 3, figsize=(24,24)) +i = 0 +for sub in Syn_diff_sc_mod: + if i < 4: + + im0 = ax[i,0].imshow(Syn_diff_sc[sub], cmap='bwr', vmax = vmax_ls[0], vmin=vmin_ls[0]) + im1 = ax[i,1].imshow(Syn_diff_sc_mod[sub], cmap='bwr', vmax = vmax_ls[1], vmin=vmin_ls[1]) + im2 = ax[i,2].imshow(Syn_diff_fc_sim[sub] - np.diag(np.diag(Syn_diff_fc_sim[sub])), cmap='bwr', vmax = vmax_ls[2], vmin=vmin_ls[2]) + + #ax[i,0].set_title(titles[i] +': initial SC', fontsize= 20) + ax[i,0].tick_params(labelsize= 15) + #ax[i,0].tick_params(axis='y', labelsize= 15) + #ax[i,1].set_title(titles[i] +': modified SC', fontsize= 20) + ax[i,1].tick_params(labelsize= 15) + #ax[i,1].tick_params(axis='y', labelsize= 15) + #ax[i,2].set_title(sub +': sim FC', fontsize= 20) + ax[i,2].tick_params(labelsize= 15)# + #ax[i,2].tick_params(axis='y', labelsize= 15) + ax[i,0].set_xticklabels([]) + ax[i,0].set_yticklabels([]) + ax[i,1].set_xticklabels([]) + ax[i,1].set_yticklabels([]) + ax[i,2].set_xticklabels([]) + ax[i,2].set_yticklabels([]) + i += 1 + else: + break +cbar0 = fig.colorbar(im0, ax = ax.T[0], location='bottom', aspect=10, fraction = 0.046, pad = 0.02) +cbar0.set_ticks(delta+np.round(np.linspace(vmin_ls[0], vmax_ls[0], 5), 3)) +cbar0.set_ticklabels(delta+ np.round(np.linspace(vmin_ls[0], vmax_ls[0], 5),3)) + +cbar1 = fig.colorbar(im1, ax = ax.T[1], location='bottom', aspect=10,fraction = 0.046, pad = 0.02) +cbar1.set_ticks(delta+np.round(np.linspace(vmin_ls[1], vmax_ls[1], 5), 3)) +cbar1.set_ticklabels(delta+np.round(np.linspace(vmin_ls[1], vmax_ls[1], 5), 3)) + +cbar2 = fig.colorbar(im2, ax = ax.T[2], location='bottom', aspect=10,fraction = 0.046, pad = 0.02) +cbar2.set_ticks(delta+np.round(np.linspace(vmin_ls[2], vmax_ls[2], 5), 3)) +cbar2.set_ticklabels(delta+np.round(np.linspace(vmin_ls[2], vmax_ls[2], 5),3)) +plt.show() +# fig.savefig(base_dir+'Syn_sc_bestfit_differentInits.png') + +fig, ax = plt.subplots(2,2, figsize=(40,28)) + +mask = np.tril_indices(66,-1) +data_dict = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] + +para_name = ['g', + 'gEE', + 'gIE', + 'gEI'] + +for sub in Syn_paras: + + data_dict['g'].append(Syn_paras[sub][0].mean()) + data_dict['gEE'].append(Syn_paras[sub][1].mean()) + data_dict['gIE'].append(Syn_paras[sub][2].mean()) + data_dict['gEI'].append(Syn_paras[sub][3].mean()) + +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) + +for i in range(4): + sns.histplot(ax=ax[i//2, i%2], data=data_dict[para_name[i]], kde=True, line_kws={'linewidth':10}) + ax[i//2, i%2].set_title('hist of ' + para_name[i] +' across subjects', fontsize='40') + ax[i//2, i%2].tick_params(labelsize=40) + ax[i//2, i%2].set_xlabel(para_name[i], fontsize=40) + +plt.show() +# fig.savefig(base_dir+'Syn_hist_modelparas.png') + +sc_m = Syn_sc_mod[sub] + +fig, ax = plt.subplots(1,3, figsize=(12,8), sharex=True, sharey=True) +ax[0].scatter(fc_true[mask], sc_m[mask], c='r') +#ax[0].set_title('The modified SC against FC in elements', fontsize=20) +ax[0].tick_params(labelsize = 20) +ax[0].set_xlabel('FC weights', fontsize = 20) +ax[0].set_ylabel('Modified SC weights', fontsize = 20) +ax[0].set_ylim(0,0.18) +ax[1].scatter(fc_true[mask], sc_true[mask]) +#ax[1].set_title('The SC against FC in elements', fontsize=20) +ax[1].set_ylim(0,0.18) +ax[1].tick_params(labelsize = 20) +ax[1].set_xlabel('FC weights', fontsize = 20) +ax[1].set_ylabel('SC weights', fontsize = 20) +ax[2].scatter(fc_true[mask], Syn_sc[sub][mask]) +#ax[2].set_title('The initial SC against FC in elements', fontsize=20) +ax[2].set_ylim(0,0.18) +ax[2].tick_params(labelsize = 20) +ax[2].set_xlabel('FC weights', fontsize = 20) +ax[2].set_ylabel('SC initial weights', fontsize = 20) +ax[0].grid() +ax[1].grid() +ax[2].grid() + +fig.tight_layout() +plt.show() + +sub = 'sub_25' + +sc = Syn_sc[sub] +sc_m = Syn_sc_mod[sub] +fig, ax = plt.subplots(1,1, figsize=(12,8)) +ax.scatter(fc_true[mask], np.corrcoef(Syn_ts_sim[sub])[mask], c='r') +#ax.set_title('The simulated FC against FC in elements',fontsize=20) +ax.tick_params(labelsize=30) +ax.set_xlabel('FC weights', fontsize=30) +ax.set_ylabel('simFC weights', fontsize=30) +ax.grid() +#ax.set_ylim(0,0.18) + +## CS add in +import seaborn as sns +import matplotlib.pyplot as plt +import numpy as np + +sub = 'sub_25' + +# Extract data +empirical_fc = fc_true[mask] # (Synthetic) empirical FC values +simulated_fc = np.corrcoef(Syn_ts_sim[sub])[mask] # Fitted FC values + +# Create histogram plot +fig, ax = plt.subplots(figsize=(12, 8)) + +# Plot empirical FC in blue +sns.histplot(empirical_fc, bins=30, color='b', label="Empirical FC", alpha=0.6) + +# Plot simulated FC in red +sns.histplot(simulated_fc, bins=30, color='r', label="Fitted FC", alpha=0.6) + +# Customize the plot +ax.tick_params(labelsize=25) +ax.set_xlabel("FC weights", fontsize=30) +ax.set_ylabel("Count", fontsize=30) +ax.legend(fontsize=20) +ax.grid() + +# Show the plot +plt.show() + +start_time = time.time() +fig, ax = plt.subplots(1,1, figsize=(12,8)) + +mask = np.tril_indices(66,-1) +data_dict = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] +"""data_dict['cA'] = [] +data_dict['cB'] = [] +data_dict['cC'] = []""" +data_dict['R'] = [] +data_dict['corr_test'] = [] +para_name = ['g', + 'gEE', + 'gIE', + 'gEI', + 'R', 'corr_test']#'cA', 'cB', 'cC', +for sub in Syn_paras: + + data_dict['g'].append(Syn_paras[sub][0].mean()) + data_dict['gEE'].append(Syn_paras[sub][1].mean()) + data_dict['gIE'].append(Syn_paras[sub][2].mean()) + data_dict['gEI'].append(Syn_paras[sub][3].mean()) + """data_dict['cA'].append(HCP_para[sub][-10:,4].mean()) + data_dict['cB'].append(HCP_para[sub][-10:,5].mean()) + data_dict['cC'].append(HCP_para[sub][-10:,6].mean())""" + data_dict['R'].append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_sim[sub][:,10:])[mask])[0,1]) + data_dict['corr_test'].append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_test[sub])[mask])[0,1]) +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) + +"""data_dict['cA'] = np.array(data_dict['cA']) +data_dict['cB'] = np.array(data_dict['cB']) +data_dict['cC'] = np.array(data_dict['cC'])""" +data_dict['R'] = np.array(data_dict['R']) + +data_dict['corr_test'] = np.array(data_dict['corr_test']) + +corr_paras = np.zeros((len(data_dict), len(data_dict))) +for i in range(len(data_dict)): + for j in range(len(data_dict)): + corr_paras[i,j] = np.corrcoef(data_dict[para_name[i]], data_dict[para_name[j]])[0,1] + +im = ax.imshow(corr_paras, cmap='bwr', vmin=-1) +ax.set_xticklabels(['0' ] + para_name, fontsize=30) +ax.set_yticklabels( ['0'] + para_name, fontsize=30) +ax.tick_params(labelsize=20) + +cbar = fig.colorbar(im, ax=ax) +cbar.ax.tick_params(labelsize=30) +cbar.set_ticks([-1.0, -0.5, 0.0, 0.5, 1.0]) +cbar.set_ticklabels([-1.0, -0.5, 0.0, 0.5, 1.0]) +plt.show() +# fig.savefig(base_dir+'Syn_corr_modelparas+fit.png') + +start_time = time.time() +fig, ax = plt.subplots(2,2, figsize=(16,12)) +plt.rcParams["axes.labelsize"] = 20 +mask = np.tril_indices(66,-1) +data_dict = {} +data_max = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] + +data_max['g'] = [] +data_max['gEE'] = [] +data_max['gIE'] = [] +data_max['gEI'] = [] + +para_name = ['g', + 'gEE', + 'gIE', + 'gEI'] + +data_chains = [] + + +for j in range(10): + data_dict['g'] = [] + data_dict['gEE'] = [] + data_dict['gIE'] = [] + data_dict['gEI'] = [] + + data_max['g'] = [] + data_max['gEE'] = [] + data_max['gIE'] = [] + data_max['gEI'] = [] + + numID=np.random.choice(np.arange(10), size= 30) + subs = ['sub_'+str(num) for num in numID] + for sub in subs: + + data_max['g'].append(max(np.abs(Syn_paras[sub][0].mean()), 100)) + data_max['gEE'].append(max(np.abs(Syn_paras[sub][1].mean()), 3.5)) + data_max['gIE'].append(max(np.abs(Syn_paras[sub][2].mean()), 0.42)) + data_max['gEI'].append(max(np.abs(Syn_paras[sub][3].mean()), 0.42)) + + for sub in subs: + + data_dict['g'].append((Syn_paras[sub][0].mean() - 100))#/max(data_max['g'])) + data_dict['gEE'].append((Syn_paras[sub][1].mean() - 3.5))#/max(data_max['gEE'])) + data_dict['gIE'].append((Syn_paras[sub][2].mean() - 0.42))#/max(data_max['gIE'])) + data_dict['gEI'].append((Syn_paras[sub][3].mean() - 0.42))#/max(data_max['gEI'])) + + data_dict['g'] = np.array(data_dict['g']) + data_dict['gEE'] = np.array(data_dict['gEE']) + data_dict['gIE'] = np.array(data_dict['gIE']) + data_dict['gEI'] = np.array(data_dict['gEI']) + data_chains.append(np.array(list(data_dict.values())).T) + + for i in range(4): + sns.kdeplot(ax=ax[i//2, i%2], data=data_dict[para_name[i]], linewidth=10) + #ax[i//2, i%2].set_title('hist of error of ' + para_name[i] +' across subjects', fontsize = 20) + #ax[i//2, i%2].set_xlim(-0.25,0.25) + ax[i//2, i%2].tick_params(labelsize = 15) + +plt.show() +# fig.savefig(base_dir + 'Syn_hist_10chains_withgains.png') + +R_stat(np.array(data_chains)) +fig, ax = plt.subplots(1, 1, figsize=(12,8)) +ax.bar(para_name, R_stat(np.array(data_chains)), color='g') +ax.plot(np.linspace(-0.5, 3.5, 50), 1.05*np.ones((50,)), c='r') +ax.set_ylabel('$\hat{R}$', fontsize=30) +ax.tick_params(labelsize=15) +plt.show() +# fig.savefig(base_dir + 'Syn_R_stat_withgains.png') + +"""# 4. Model fitting and training""" + +out_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP' +base_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP' +subs =sorted([sc_file[-10:-4] for sc_file in os.listdir(base_dir) if sc_file[:8] == 'weights_']) + +"""node_size = 83 +step_size = 0.05 +Tr = 2.5 + +g = 82 +gEE = 2.5 +gIE = 0.4473 +gEI = 0.501 + +sc_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/weights.txt' #CS +sc = np.loadtxt(sc_file) +sc = (sc + sc.T)*0.5 +SC = sc - np.diag(np.diag(sc)) +#SC[SC 0].mean() +1.1*SC[SC>0].std()] =0 # for connectivity 96 only +wo = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) + +test_sim = WWD_test(g, gEE, gIE, gEI, wo, step_size, node_size, Tr) + +bold_test = test_sim.generate_bold(4800) #could be any number you want """ + +# Define the file path where you want to save the result +ts_file = "bold_test_output.txt" # Change the name if needed + +# Save the array to a text file +np.savetxt(ts_file, bold_test) + +print(f"Saved bold_test results to {ts_file}") + # CS add in + +base_dir = '../wwd-model-fitting/data/Hagmann/' #cs +#ts_file = base_dir+'syn_ts.txt' CS + +#SC = np.loadtxt(SC_file) +ts = np.loadtxt(ts_file) +print(ts.shape) +sc_file = base_dir + 'weights.txt' + +start_time = time.time() + +sub = ['150423', '257542', '154936'] + +for i in range(40): + + + base_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/' + + node_size = 83 + mask = np.tril_indices(node_size, -1) + num_epoches = 20 + batch_size = 20 + step_size = 0.05 + input_size = 2 + tr = 0.75 + sub=subs[i] + print(i, sub) + sc_file = base_dir +'weights_'+sub+'.txt' + ts_file = base_dir +sub+'_rfMRI_REST1_LR_hpc200_clean__l2k8_sc33_ts.pkl'#out_dir+'sub_'+sub+'simBOLD_idt.txt'# + + if os.path.isfile(sc_file) and os.path.isfile(ts_file): + sc = np.loadtxt(sc_file) + SC =(sc+sc.T)*0.5 + + sc = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) + + + + + + ts_pd = pd.read_pickle(ts_file) + ts = ts_pd.values + #ts = np.loadtxt(ts_file) + ts =ts/np.max(ts) + fc_emp = np.corrcoef(ts.T) + # Get the WWD model module for forward in a batch. + + + model = RNNWWD(input_size, node_size, batch_size, step_size, tr, sc, True, g_mean_ini=80, g_std_ini = .1, gEE_mean_ini=2.5, gEE_std_ini = .1) + + + + # call model fit method + F = Model_fitting(model, ts, num_epoches) + + # fit data(train) + + output_train = F.train() + + + output_test = F.test(120) + plot_fit_parameters(output_train) + plot_sim_states_outputs(ts, output_test) + + sc_mod = np.zeros_like(sc) + sc_mod[mask] = output_train['gains'][-100:].mean(0) + sc_mod = sc_mod+sc_mod.T + w = (1 + np.tanh(sc_mod))*sc + w_n = 0.5*(w + w.T)/np.linalg.norm(0.5*(w + w.T)) + np.savetxt(out_dir + 'bold_test_'+ sub +'.txt', output_test['simBOLD']) + np.savetxt(out_dir + 'bold_train_'+ sub +'.txt', output_train['simBOLD']) + np.savetxt(out_dir + 'sc_mod_'+ sub +'.txt', sc_mod) + #np.savetxt(out_dir + 'sc_'+ sub +'.txt', sc) + g= output_train['g'][-100:].mean() + gEE = output_train['gEE'][-100:].mean() + gEI = output_train['gEI'][-100:].mean() + gIE = output_train['gIE'][-100:].mean() + + + np.savetxt(out_dir + 'parameters_'+ sub +'.txt', np.array([g,gEE, gIE, gEI]).T) +end_time = time.time() +print('running time is {0} \'s'.format(end_time - start_time )) + +start_time = time.time() + + +for i in range(0, 40): + sub=subs[i] + print(i, sub) + + + node_size = 83 + num_epoches = 20 + batch_size =20 + step_size = 0.05 + input_size = 2 + tr = 0.75 + sc_file = base_dir +'weights_'+sub+'.txt' + ts_file = out_dir + 'bold_test_'+ sub +'.txt' + + if os.path.isfile(sc_file) and os.path.isfile(ts_file): + + sc = np.loadtxt(sc_file) + SC =(sc+sc.T)*0.5 + + sc = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) + + + + sc_mod = np.zeros_like(sc) + + + ts = np.loadtxt(ts_file) + ts = (ts/np.max(ts))[:,:1200] + fc_emp = np.corrcoef(ts[:,:1200]) + print(fc_emp.shape) + # Get the WWD model module for forward in a batch. + + + model = RNNWWD(input_size, node_size, batch_size, step_size, tr, sc, True, g_mean_ini=100, g_std_ini = .1, gEE_mean_ini=2.5, gEE_std_ini = .1) + + + + # call model fit method + F = Model_fitting(model, ts.T, num_epoches) + + # fit data(train) + + output_train = F.train() + + + output_test = F.test(120) + plot_fit_parameters(output_train) + plot_sim_states_outputs(ts.T, output_test) + sc_mod[mask] = output_train['gains'][-100:].mean(0) + sc_mod = sc_mod+sc_mod.T + w = (1 + np.tanh(sc_mod))*sc + w_n = 0.5*(w + w.T)/np.linalg.norm(0.5*(w + w.T)) + np.savetxt(out_dir + 'bold_test_idt_'+ sub +'.txt', output_test['simBOLD']) + np.savetxt(out_dir + 'bold_train_idt_'+ sub +'.txt', output_train['simBOLD']) + np.savetxt(out_dir + 'sc_mod_idt_'+ sub +'.txt', sc_mod) + #np.savetxt(out_dir + 'sc_'+ sub +'.txt', sc) + g= output_train['g'][-100:].mean() + gEE = output_train['gEE'][-100:].mean() + gEI = output_train['gEI'][-100:].mean() + gIE = output_train['gIE'][-100:].mean() + + + np.savetxt(out_dir + 'parameters_idt_'+ sub +'.txt', np.array([g,gEE, gIE, gEI]).T) +end_time = time.time() +print('running time is {0} \'s'.format(end_time - start_time )) + +np.loadtxt(ts_file).shape + +output_test['simBOLD'].shape + +"""HCP_ts_sim ={} +HCP_ts_test ={} +HCP_paras = {} +HCP_ts = {} +HCP_sc = {} +HCP_sc_mod = {} +for i in range(40): + sub = subs[i] + + """ts_file = base_dir +sub+'_rfMRI_REST1_LR_hpc200_clean__l2k8_sc33_ts.pkl'#out_dir+'sub_'+sub+'simBOLD_idt.txt'# + + ts_pd = pd.read_pickle(ts_file) + ts = ts_pd.values + HCP_ts[sub] = ts""" + HCP_ts_test[sub] = np.loadtxt(out_dir + 'bold_test_'+ sub +'.txt') + HCP_ts_sim[sub] = np.loadtxt(out_dir + 'bold_train_'+ sub +'.txt') + + HCP_sc_mod[sub] = np.loadtxt(out_dir + 'sc_mod_'+ sub +'.txt') + + HCP_paras[sub] = np.loadtxt(out_dir + 'parameters_'+ sub +'.txt') + +np.save(out_dir +'HCP_ts_sim.npy', HCP_ts_sim) +np.save(out_dir + 'HCP_ts_test.npy', HCP_ts_test) +#np.save(out_dir + 'HCP_ts.npy', HCP_ts) +np.save(out_dir + 'HCP_sc_mod.npy', HCP_sc_mod) +np.save(out_dir + 'HCP_fitparas.npy', HCP_paras)""" + +HCP_ts_sim ={} +HCP_ts_test ={} +HCP_paras = {} + +HCP_sc_mod = {} +for i in range(40): + sub = subs[i] + + + HCP_ts_test[sub] = np.loadtxt(out_dir + 'bold_test_idt_'+ sub +'.txt') + HCP_ts_sim[sub] = np.loadtxt(out_dir + 'bold_train_idt_'+ sub +'.txt') + + HCP_sc_mod[sub] = np.loadtxt(out_dir + 'sc_mod_idt_'+ sub +'.txt') + + HCP_paras[sub] = np.loadtxt(out_dir + 'parameters_idt_'+ sub +'.txt') + +for i in range(3): + sub = subs_s[i] + + + HCP_ts_test[sub] = np.loadtxt(out_dir + 'bold_test_idt_'+ sub +'.txt') + HCP_ts_sim[sub] = np.loadtxt(out_dir + 'bold_train_idt_'+ sub +'.txt') + + HCP_sc_mod[sub] = np.loadtxt(out_dir + 'sc_mod_idt_'+ sub +'.txt') + + HCP_paras[sub] = np.loadtxt(out_dir + 'parameters_idt_'+ sub +'.txt') +np.save(out_dir +'HCP_ts_sim_idt.npy', HCP_ts_sim) +np.save(out_dir + 'HCP_ts_test_idt.npy', HCP_ts_test) +np.save(out_dir + 'HCP_sc_mod_idt.npy', HCP_sc_mod) + + +np.save(out_dir + 'HCP_fitparas_idt.npy', HCP_paras) + +"""### Analysis of the HCP data model-fitting +Approximate run duration: 10-20 seconds + +The model-fitting methods on the data from the Human Connectome Project (HCP) let us view how well the parameters fit the model. + +The HCP (a.k.a."original" or main HCP, HCP Young Adult, HCP-YA) maps the healthy human connectome by collecting and freely distributing neuroimaging and behavioral data on 1,200 normal young adults, aged 22-35. This notebook contains the post-analysis in which we determine how well the parameters fit the model using the model-fitting method. + +After fitting the model to the empirical BOLD signal, one may get a set of model parameters that fit the model well. We can use the forward model and the fitted model parameters to generate simulated BOLD data, fit the model again to the simulated BOLD signal with the fitted model parameters (true values), and the new model-fitted data are saved at the files with the `idt` suffix. + +In this analysis, we demonstrate how well they fit. + +#### Initialize variables + +Import the HCP data and corresponding connectivity matrix files. +""" + +#base_dir = '../data/HCP/Outputs/' # base directory for test HCP data +base_dir = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/' # base directory for test HCP data +fc_file = base_dir + 'HCP_ts.npy' # HCP functional connectivity data from the HCP dataset +sc_file = base_dir + 'HCP_sc.npy' # HCP structural connectivity data +fc_sim_file = base_dir + 'HCP_ts_sim.npy' # HCP functional connectivity data from the simulations +fc_sim_idt_file = base_dir + 'HCP_ts_sim_idt.npy' # HCP simulated functional connectivity file with the + # model-identified parameters +sc_mod_file = base_dir + 'HCP_sc_mod.npy' # HCP modified structural connectivity data +sc_mod_idt_file = base_dir + 'HCP_sc_mod_idt.npy' # newly model-fitted structural connectivity data +paras_file = base_dir + 'HCP_fitparas.npy' # fitted parameters file +paras_idt_file = base_dir + 'HCP_fitparas_idt.npy'# new model-fitted data +fc_test_file = base_dir + 'HCP_ts_test.npy' # test data from HCP +fc_test_idt_file = base_dir + 'HCP_ts_test_idt.npy' # test data from HCP + +HCP_fc_data = np.load(fc_file, allow_pickle=True) +HCP_sc_data = np.load(sc_file, allow_pickle=True) +HCP_fc_sim_data = np.load(fc_sim_file, allow_pickle=True) +HCP_fc_sim_idt_data = np.load(fc_sim_idt_file, allow_pickle=True) +HCP_sc_mod_data = np.load(sc_mod_file, allow_pickle=True) +HCP_sc_mod_idt_data = np.load(sc_mod_idt_file, allow_pickle=True) +HCP_para_data = np.load(paras_file, allow_pickle=True) +HCP_para_idt_data = np.load(paras_idt_file, allow_pickle=True) +HCP_fc_sim_data = np.load(fc_sim_file, allow_pickle=True) +HCP_fc_test_data = np.load(fc_test_file, allow_pickle=True) +HCP_fc_test_idt_data = np.load(fc_test_idt_file, allow_pickle=True) + +HCP_ts = HCP_fc_data.item() +HCP_sc = HCP_sc_data.item() +HCP_ts_sim = HCP_fc_sim_data.item() +HCP_ts_sim_idt = HCP_fc_sim_idt_data.item() +HCP_sc_mod = HCP_sc_mod_data.item() +HCP_sc_mod_idt = HCP_sc_mod_idt_data.item() +HCP_para = HCP_para_data.item() +HCP_para_idt = HCP_para_idt_data.item() +HCP_ts_test = HCP_fc_test_data.item() +HCP_ts_test_idt = HCP_fc_test_idt_data.item() + +i =0 +for sub in HCP_para: + if i < 43: + print(sub) + fig, ax = plt.subplots(1, 4, figsize=(12,8)) + im0 = ax[0].imshow(HCP_sc[sub], cmap='bwr') + sc_mod = HCP_sc_mod[sub] + sc = HCP_sc[sub] + w = (1 + np.tanh(sc_mod))*sc + w_n = 0.5*(w + w.T)/np.linalg.norm(0.5*(w + w.T)) + im1 = ax[1].imshow(w_n, cmap='bwr') + im2 = ax[2].imshow(np.corrcoef((HCP_ts[sub]).T)-np.diag(np.diag(np.corrcoef((HCP_ts[sub]).T))), cmap='bwr', vmax =1) #-HCP_ts[sub].mean(1))- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T-HCP_ts[sub].mean(1)))) + im3 = ax[3].imshow(np.corrcoef(HCP_ts_test[sub])-np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]))), cmap='bwr',vmax = 1) #-HCP_ts_test[sub].mean(0))- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]-HCP_ts_test[sub].mean(0)))) + plt.show() + i +=1 + +"""#### Plot the connectivity weight matrices""" + +fig, ax = plt.subplots(3, 4, figsize=(20,16)) + +vmax_ls = [0,-50,0, 0] +vmin_ls = [1,50,1, 1] +delta= np.array([0,0,0,0,-0.001]) +i = 0 +#for sub in ['562446', '257542', '154936']: +for sub in ['150423', '257542', '154936']: # CS choose + if i < 4: + if HCP_sc[sub].max() > vmax_ls[0]: + vmax_ls[0] = HCP_sc[sub].max() + if HCP_sc_mod[sub].max() > vmax_ls[1]: + + vmax_ls[1] = HCP_sc_mod[sub].max() + if (np.corrcoef(HCP_ts[sub].T)- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T)))).max() > vmax_ls[2]: + vmax_ls[2] = (np.corrcoef(HCP_ts[sub].T)- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T)))).max() + if (np.corrcoef(HCP_ts_test[sub])- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub])))).max() > vmax_ls[3]: + vmax_ls[3] = (np.corrcoef(HCP_ts_test[sub])- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub])))).max() + if HCP_sc[sub].min() < vmin_ls[0]: + vmin_ls[0] = HCP_sc[sub].min() + if HCP_sc_mod[sub].min() < vmin_ls[1]: + + vmin_ls[1] = HCP_sc_mod[sub].min() + if (np.corrcoef(HCP_ts[sub].T)- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T)))).min() < vmin_ls[2]: + vmin_ls[2] = (np.corrcoef(HCP_ts[sub].T)- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T)))).min() + if(np.corrcoef(HCP_ts_test[sub])- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub])))).min() < vmin_ls[3]: + vmin_ls[3] = (np.corrcoef(HCP_ts_test[sub])- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub])))).min() + i += 1 +i_st = 0 +i = 0 +delta = np.array([0.001,0,0,0,0.0]) +for sub in ['150423', '257542', '154936']: # CS choose +#for sub in ['562446', '257542', '154936'] : + if (i > i_st-1 ) and (i < i_st + 3): + im0 = ax[i-i_st,0].imshow(HCP_sc[sub], cmap='bwr') + sc_mod = HCP_sc_mod[sub] + sc = HCP_sc[sub] + + im1 = ax[i-i_st,1].imshow(sc_mod, cmap='bwr') + im2 = ax[i-i_st,2].imshow(np.corrcoef((HCP_ts[sub]).T)-np.diag(np.diag(np.corrcoef((HCP_ts[sub]).T))), cmap='bwr', vmax = 0.8, vmin=-0.2) #-HCP_ts_test[sub].mean(0))- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]-HCP_ts_test[sub].mean(0)))) + im3 = ax[i-i_st,3].imshow(np.corrcoef(HCP_ts_test[sub])-np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]))), cmap='bwr', vmax= 0.8, vmin=-0.2)#-HCP_ts_test[sub].mean(0))- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]-HCP_ts_test[sub].mean(0))) + #ax[i-i_st,0].set_title('s'+sub +': empirical SC', fontsize=20) + ax[i-i_st,0].set_xticklabels([]) + ax[i-i_st,0].set_yticklabels([]) + #ax[i-i_st,1].set_title('s'+sub +': modified SC', fontsize=20) + ax[i-i_st,1].set_xticklabels([]) + ax[i-i_st,1].set_yticklabels([]) + #ax[i,0].tick_params(axis='y', labelsize=15) + #ax[i-i_st,2].set_title('s'+sub +': empirical FC', fontsize=20) + ax[i-i_st,2].set_xticklabels([]) + ax[i-i_st,2].set_yticklabels([]) + #ax[i,1].tick_params(axis='y', labelsize=15) + #ax[i-i_st,3].set_title('s'+sub +' simulated FC', fontsize=20) + ax[i-i_st,3].set_xticklabels([]) + ax[i-i_st,3].set_yticklabels([]) + #ax[i,2].tick_params(axis='y', labelsize= 15) + i += 1 + elif i < i_st: + i += 1 + else: + break + +cbar0 = fig.colorbar(im0, + ax=ax.T[0], + location='bottom', + aspect=10, + fraction=0.046, + pad=0.02) + +cbar1 = fig.colorbar(im1, + ax=ax.T[1], + location='bottom', + aspect=10, + fraction=0.046, + pad=0.02) + +cbar2 = fig.colorbar(im2, + ax=ax.T[2], + location='bottom', + aspect=10, + fraction=0.046, + pad=0.02) + +cbar3 = fig.colorbar(im3, + ax=ax.T[3], + location='bottom', + aspect=10, + fraction=0.046, + pad=0.02) + +cbar0.set_ticks(0.001*np.floor(1000*np.linspace(vmin_ls[0], vmax_ls[0], 5))) + +cbar1.set_ticks(0.001*np.floor(1000*np.linspace(vmin_ls[1], vmax_ls[1], 5))) +cbar2.set_ticks(0.001*np.floor(1000*np.linspace(-0.2, 0.8, 5))) +cbar3.set_ticks(0.001*np.floor(1000*np.linspace(-0.2, 0.8, 5))) +plt.show() +# fig.savefig(base_dir+'HCP_sc_fc_bestfit_3subjects.png') + +"""#### Model parameter fitting estimation""" + +start_time = time.time() +mask = np.tril_indices(83,-1) + +data_dict = {} +data_max = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] + +data_max['g'] = [] +data_max['gEE'] = [] +data_max['gIE'] = [] +data_max['gEI'] = [] + +para_name = ['g', 'gEE', 'gIE', 'gEI'] + +""" +Across the HCP test subjects, take the maximum of the absolute value of the mean of the parameter +values identified for the model and the same for the parameter values. Append the two values and +add them to the `data_max` dictionary for the coupling strength parameters. +""" + +for sub in HCP_para: + + data_max['g'].append(max(np.abs(HCP_para_idt[sub][0].mean()), np.abs(HCP_para[sub][0].mean()))) + data_max['gEE'].append(max(np.abs(HCP_para_idt[sub][1].mean()), np.abs(HCP_para[sub][1].mean()))) + data_max['gIE'].append(max(np.abs(HCP_para_idt[sub][2].mean()), np.abs(HCP_para[sub][2].mean()))) + data_max['gEI'].append(max(np.abs(HCP_para_idt[sub][3].mean()), np.abs(HCP_para[sub][3].mean()))) + +""" +Take these parameter values and use these maximum values to normalize the range of parameter values. +""" + +for sub in HCP_para: + + data_dict['g'].append((HCP_para_idt[sub][0].mean() - HCP_para[sub][0].mean())/max(data_max['g'])) + data_dict['gEE'].append((HCP_para_idt[sub][1].mean() - HCP_para[sub][1].mean())/max(data_max['gEE'])) + data_dict['gIE'].append((HCP_para_idt[sub][2].mean() - HCP_para[sub][2].mean())/max(data_max['gIE'])) + data_dict['gEI'].append((HCP_para_idt[sub][3].mean() - HCP_para[sub][3].mean())/max(data_max['gEI'])) + +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) + +fig, ax = plt.subplots(2,2, figsize=(12,8), sharex=True) +plt.rcParams["axes.labelsize"] = 40 +for i in range(4): + sns.kdeplot(ax=ax[i//2, i%2], data=data_dict[para_name[i]], color='r', linewidth=10) + #ax[i//2, i%2].set_title('hist of error of ' + para_name[i] +' across subjects', fontsize = 20) + ax[i//2, i%2].set_xlim(-0.5,0.5) + ax[i//2, i%2].tick_params(labelsize = 30) + ax[i//2, i%2].grid() + if i%2 == 1: + ax[i//2, i%2].set_ylabel('') +fig.tight_layout() +plt.show() +# fig.savefig(base_dir+'HCP_hist_modelparaserror_normalized_pt.png') + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 25 +sns.kdeplot(data=pd.DataFrame(data_dict), linewidth=4) +#ax.set_title('hist of error of parameters across subjects', fontsize =25) +ax.tick_params(labelsize=20) +ax.set_xlim(-0.2,0.2) +ax.set_ylim(0,50) +ax.grid() +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize='50') +plt.show() +# fig.savefig(base_dir+'HCP_hist_modelparaserror_normalized_in1fig_pt.png') + +fig, ax = plt.subplots(2,2, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 40 +mask = np.tril_indices(83,-1) +data_dict = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] +data_dict['cA'] = [] +data_dict['cB'] = [] +data_dict['cC'] = [] +para_name = ['g', 'gEE', 'gIE', 'gEI'] #, 'cA', 'cB', 'cC'] + +""" +Get the mean of the parameters in the dictionary of data. +""" + +for sub in HCP_para: + data_dict['g'].append(HCP_para[sub][0].mean()) + data_dict['gEE'].append(HCP_para[sub][1].mean()) + data_dict['gIE'].append(HCP_para[sub][2].mean()) + data_dict['gEI'].append(HCP_para[sub][3].mean()) + """ + data_dict['cA'].append(HCP_para[sub][-10:,4].mean()) + data_dict['cB'].append(HCP_para[sub][-10:,5].mean()) + data_dict['cC'].append(HCP_para[sub][-10:,6].mean()) + """ +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) + +""" +data_dict['cA'] = np.array(data_dict['cA']) +data_dict['cB'] = np.array(data_dict['cB']) +data_dict['cC'] = np.array(data_dict['cC']) +""" +for i in range(4): + sns.kdeplot(ax =ax[i//2, i%2], data= data_dict[para_name[i]], color='r', linewidth = 10) + #ax[i//2, i%2].set_title('hist of ' + para_name[i] +' across subjects', fontsize = 20) + ax[i//2, i%2].tick_params(labelsize = 20) + ax[i//2, i%2].set_xlabel(para_name[i], fontsize = 20) + ax[i//2, i%2].grid() + if i%2 == 1: + ax[i//2, i%2].set_ylabel('') +fig.tight_layout() +plt.show() +#fig.savefig(base_dir+'HCP_hist_modelparas.png') + +fig, ax = plt.subplots(2,2, figsize=(12,8), sharex=True, sharey=True) +plt.rcParams["axes.labelsize"] = 20 +i = 0 +mask = np.tril_indices(83,-1) +for sub in ['198653', '257542', '154936', '280941']: + if i < 4: + data_dict = {} + data_dict['empFC'] = np.corrcoef(HCP_ts[sub].T)[mask] + data_dict['simFC'] = np.corrcoef(HCP_ts_sim[sub][:,10:])[mask] + data_dict['testFC'] = np.corrcoef(HCP_ts_test[sub])[mask] + sns.kdeplot(ax =ax[i//2, i%2], data=pd.DataFrame(data_dict), linewidth=10) + #ax[i//2, i%2].set_title('s'+sub, fontsize= 30) + ax[i//2, i%2].tick_params(labelsize=20) + ax[i//2, i%2].grid() + plt.setp(ax[i//2, i%2].get_legend().get_texts(), fontsize='25') + plt.setp(ax[i//2, i%2].get_legend().get_title(), fontsize ='40') + if i%2 == 1: + ax[i//2, i%2].set_ylabel('') + i += 1 +fig.tight_layout() +plt.show() +# fig.savefig(base_dir+'HCP_hist_FCs_pt.png') + +fig, ax = plt.subplots(2,2,figsize=(12,8), sharex=True, sharey=True) +plt.rcParams["axes.labelsize"] = 20 +i = 0 +mask = np.tril_indices(83, -1) +for sub in ['198653', '257542', '154936', '280941']: + """ + The difference between original and modified structural connectivity weights. + """ + if i < 4: + data_dict = {} + data_dict['originalSC'] = HCP_sc[sub][mask] + sc_mod = HCP_sc_mod[sub] + sc = HCP_sc[sub] + w = (1 + np.tanh(sc_mod))*sc + w_n = 0.5*(w + w.T)/np.linalg.norm(0.5*(w + w.T)) + data_dict['modifiedSC'] = sc_mod[mask] + sns.kdeplot(ax=ax[i//2, i%2], data=pd.DataFrame(data_dict), linewidth=10) + #ax[i//2, i%2].set_title('s'+sub, fontsize = 30) + ax[i//2, i%2].tick_params(labelsize = 25) + ax[i//2, i%2].grid() + plt.setp(ax[i//2, i%2].get_legend().get_texts(), fontsize='25') + plt.setp(ax[i//2, i%2].get_legend().get_title(), fontsize='40') + + if i%2 == 1: + ax[i//2, i%2].set_ylabel('') + i += 1 +fig.tight_layout() +plt.show() +# fig.savefig(base_dir+'HCP_hist_SCs_pt.png') + +corr_sim = {} +corr_sim['test'] = [] +corr_sim['train'] = [] +mask = np.tril_indices(83,-1) +for sub in HCP_ts: + corr_sim['train'].append(np.corrcoef(np.corrcoef(HCP_ts_sim[sub][:,10:])[mask], np.corrcoef(HCP_ts[sub].T)[mask])[0,1]) + corr_sim['test'].append(np.corrcoef(np.corrcoef(HCP_ts_test[sub])[mask], np.corrcoef(HCP_ts[sub].T)[mask])[0,1]) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 40 +sns.kdeplot(data=corr_sim, linewidth=10) +#ax.set_title('hist of Pearson correlation between the fitting and empirical BOLD', fontsize=20) +ax.tick_params(labelsize=25) +ax.grid() +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize='50') + +# plt.savefig(base_dir+'HCP_hist_fittingcorrelation_pt.png') + +# cross valid for identifiability... trained on the fist 2400 and tested on the last 2400 data points +corr_sim = {} +corr_sim['test'] = [] +corr_sim['train'] = [] +mask = np.tril_indices(83,-1) +for sub in HCP_ts_test: + corr_sim['train'].append(np.corrcoef(np.corrcoef(HCP_ts_sim_idt[sub][:,10:])[mask], np.corrcoef(HCP_ts_test[sub][:,:1200])[mask])[0,1]) + corr_sim['test'].append(np.corrcoef(np.corrcoef(HCP_ts_test_idt[sub])[mask], np.corrcoef(HCP_ts_test[sub][:,1200:])[mask])[0,1]) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 40 +sns.kdeplot(data=corr_sim, linewidth=10) +#ax.set_title('hist of Pearson correlation between the fitting and empirical BOLD', fontsize=20) +ax.tick_params(labelsize=25) +ax.grid() +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize='50') + +fig, ax = plt.subplots(1,1, figsize=(12,8)) + +mask = np.tril_indices(83,-1) +data_dict = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] +""" +data_dict['cA'] = [] +data_dict['cB'] = [] +data_dict['cC'] = [] +""" +data_dict['$\mathregular{R_{train}}$'] = [] +data_dict['$\mathregular{R_{test}}$'] = [] +para_name = ['g', + 'gEE', + 'gIE', + 'gEI', + '$\mathregular{R_{train}}$', + '$\mathregular{R_{test}}$'] #'cA', 'cB', 'cC', + +for sub in HCP_ts: + data_dict['g'].append(HCP_para[sub][0].mean()) + data_dict['gEE'].append(HCP_para[sub][1].mean()) + data_dict['gIE'].append(HCP_para[sub][2].mean()) + data_dict['gEI'].append(HCP_para[sub][3].mean()) + """ + data_dict['cA'].append(HCP_para[sub][-10:,4].mean()) + data_dict['cB'].append(HCP_para[sub][-10:,5].mean()) + data_dict['cC'].append(HCP_para[sub][-10:,6].mean()) + """ + data_dict['$\mathregular{R_{train}}$'].append(np.corrcoef(np.corrcoef(HCP_ts[sub].T)[mask], np.corrcoef(HCP_ts_sim[sub][:,10:])[mask])[0,1]) + data_dict['$\mathregular{R_{test}}$'].append(np.corrcoef(np.corrcoef(HCP_ts[sub].T)[mask], np.corrcoef(HCP_ts_test[sub])[mask])[0,1]) +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) +""" +data_dict['cA'] = np.array(data_dict['cA']) +data_dict['cB'] = np.array(data_dict['cB']) +data_dict['cC'] = np.array(data_dict['cC']) +""" +data_dict['$\mathregular{R_{train}}$'] = np.array(data_dict['$\mathregular{R_{train}}$']) +data_dict['$\mathregular{R_{test}}$'] = np.array(data_dict['$\mathregular{R_{test}}$']) + +corr_paras = np.zeros((len(data_dict), len(data_dict))) +for i in range(len(data_dict)): + for j in range(len(data_dict)): + corr_paras[i,j] = np.corrcoef(data_dict[para_name[i]], data_dict[para_name[j]])[0,1] + +im = ax.imshow(corr_paras, cmap='bwr') +ax.set_xticklabels(['0' ] + para_name, fontsize=20) +ax.set_yticklabels(['0'] + para_name, fontsize=20) +cbar = fig.colorbar(im, ax=ax) +cbar.ax.tick_params(labelsize=30) +cbar.set_ticks(np.linspace(-1,1,5)) +cbar.set_ticklabels(np.linspace(-1,1,5)) +plt.show() +# fig.savefig(base_dir + 'HCP_corr_modelparas+fit_pt.png') + +"""#### Load SC data""" + +"""l2k8_to_aparc_sort_idx = np.array([22, 21, 20, 23, 18, 17, 19, 14, 24, 30, 16, 28, 10, 25, 15, 13, 29, + 32, 31, 9, 26, 33, 27, 8, 5, 12, 7, 4, 3, 11, 0, 1, 6, 2, 34, 35, + 64, 62, 63, 65, 60, 59, 61, 56, 66, 72, 58, 70, 67, 71, 52, 57, 55, + 74, 73, 51, 68, 75, 69, 50, 47, 54, 49, 46, 53, 42, 43, 45, 48, 44, 76,77])""" + +l2k8_to_aparc_sort_idx = np.array([82, 30, 12, 8, 82, 20, 26, 24, 18, 28, 14, 22, 0, 23, 3, 29, 25, + 10, 5, 1, 4, 21, 15, 13, 9, 19, 11, 6, 7, 17, 31, 16, 2, 27, + 32, 33, 82, 71, 53, 49, 82, 61, 67, 65, 59, 69, 55, 63, 41, 64, 44, + 70, 66, 51, 46, 42, 45, 62, 56, 54, 50, 60, 52, 47, 48, 58, 72, 57, + 43, 68, 73, 74]) + +# Load the data from the .npy files. +l2k8_w = np.loadtxt('/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/hcp_grpavg_l2k8_sc33_weights.npy') # structural connectivity weights +l2k8_tl = np.loadtxt('/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/hcp_grpavg_l2k8_sc33_tractlengths.npy') # tract lengths between regions + +# Convert the files to pandas DataFrames. +df_l2k8_w = pd.DataFrame(l2k8_w) +df_l2k8_tl = pd.DataFrame(l2k8_tl) + +df_l2k8_w *= ((np.eye(83)*-1) + 1) # Multiply by a matrix of all 1s with 0s on the diagonal. +df_l2k8_wlog = np.log1p(df_l2k8_w) # Obtain log(1+x) for the weights x for the 83 x 83 matrix of regions. + +df_l2k8_wlog_divmax = df_l2k8_wlog.copy() # Copy the DataFrame. +df_l2k8_wlog_divmax /= df_l2k8_wlog.max() # Divide the DataFrame by its maximum amount to normalize it. + +_df_w = df_l2k8_wlog_divmax.copy() +_df_w[_df_w<0.2] = 0 + +# Use these regions of the brain as focal nodes in representing the DMN. + +l2k8_lV1_idx = 21 # pericalcarine +l2k8_lM1_idx = 24 # precentral +l2k8_lA1_idx = 30 # superior temporal +l2k8_lS1_idx = 22 # postcentral +l2k8_lDLPFC_idx = 12 # lateral orbitofrontal +l2k8_lVLPFC_idx = 19 # pars orbitalis + +focal_nodes = [l2k8_lV1_idx, + l2k8_lA1_idx, + l2k8_lM1_idx, + l2k8_lS1_idx, + l2k8_lDLPFC_idx, + l2k8_lVLPFC_idx] + +focal_node_names = ['lV1', + 'lA1', + 'lM1', + 'lS1', + 'lDLPFC', + 'lVLPFC'] + +# Various surface representations for left and right hemisphere +# With the inflated surface representations, you can fully see the sulci. +lhi_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/fsav5_lh.inflated' +rhi_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/fsav5_rh.inflated' + +# The pial surface is used to calculate cortical gray matter volume and +# should accurately follow the boundary between the gray matter and the Cerebrospinal fluid (CSF). +lhp_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/fsav5_lh.pial' +rhp_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/fsav5_rh.pial' + +# Annotatation files for the parcellated regions +lh_annot_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/fsav5_lh.aparc.annot' +rh_annot_file = '/content/drive/MyDrive/wwd-model-fitting-main/data/HCP/fsav5_rh.aparc.annot' + +# Return the vertex coordinates (rhi_vtx) and mesh triangles (rhi_tri) for the files +lhi_vtx, lhi_tri = nib.freesurfer.read_geometry(lhi_file) +rhi_vtx, rhi_tri = nib.freesurfer.read_geometry(rhi_file) + +lhp_vtx, lhp_tri = nib.freesurfer.read_geometry(lhp_file) +rhp_vtx, rhp_tri = nib.freesurfer.read_geometry(rhp_file) + +lh_annot = nib.freesurfer.read_annot(lh_annot_file) +rh_annot = nib.freesurfer.read_annot(rh_annot_file) + +vtx_l2k8_lh = lhp_vtx +tri_l2k8_lh = lhp_tri +#vtx_l2k8_lh = lhi_vtx +#tri_l2k8_lh = lhi_tri +rm_l2k8_lh = lh_annot[0] + +vtx_l2k8_rh = rhp_vtx +tri_l2k8_rh = rhp_tri +#vtx_l2k8_rh = rhi_vtx +#tri_l2k8_rh = rhi_tri +rm_l2k8_rh = rh_annot[0] + +# Concatenate them together for the regions. +vtx_l2k8_lhrh = np.concatenate([vtx_l2k8_lh,vtx_l2k8_rh], axis=0) +tri_l2k8_lhrh = np.concatenate([tri_l2k8_lh,tri_l2k8_rh+10242], axis=0) + +rm_l2k8_rh_p1 = rm_l2k8_rh + 36 +rm_l2k8_rh_p1[rm_l2k8_rh == 0] = 0 + +rm_l2k8_lhrh = np.concatenate([rm_l2k8_lh, rm_l2k8_rh_p1], axis=0) +rm_l2k8_lhrh + +hemi = np.ones_like(rm_l2k8_lhrh) +hemi[:10242] = 0 +hemi + + + +"""## Plot + +Generate the plots of the default mode network structure. +""" + +# functional connectivity values across regions for three different subjects +# The number coreresponds to the different columns of the FC (e.g., fc_54 has the 54th column). +fc_54 = np.array([0.07682509, -0.02128439, -0.03384062, 0.05320176, 0.05471503, + 0.13934857, 0.20621781, 0.17206639, 0.10512446, 0.10776393, + 0.29780534, 0.14799411, 0.30186789, 0.43196902, 0.33156828, + 0.25412462, 0.30922221, 0.12608491, 0.27343265, 0.41269201, + 0.25245327, 0.09809735, -0.11438562, 0.20858446, 0.07015369, + 0.08555236, 0.08659411, 0.04477133, 0.08317634, 0.11445271, + 0.12312393, 0.16196342, 0.06828829, 0.26020704, 0.12544554, + 0.05920674, 0.03225544, 0.01331936, -0.0215654 , 0.07715058, + -0.0389014 , 0.17158385, 0.0094781 , 0.02607858, 0.11923761, + 0.06269088, 0.10788386, 0.36012597, 0.22229945, 0.21825426, + 0.02970901, 0.24230718, 0.15896073, 0.40808317, 1. , + 0.33747972, 0.23850011, 0.30159616, 0.1130841 , 0.2969263 , + 0.47867023, 0.19959395, 0.06715806, -0.09413906, 0.20176679, + 0.10768835, 0.02944969, 0.04499624, -0.03408385, 0.17570769, + 0.16808477, 0.06259806, 0.26563371, 0.14068841, 0.27621275, + 0.1306482 , 0.09807702, 0.09850341, 0.05073915, 0.03550151, + 0.17014859, 0.06636479, 0.17870341]) + +fc_60 = np.array([0.15569094, -0.03885425, 0.04890683, 0.16171774, 0.04725465, + 0.07753457, 0.30875898, 0.35187976, 0.34662134, 0.09912189, + 0.29338462, 0.25739101, 0.17445902, 0.4670589 , 0.62694002, + 0.25998084, 0.25925487, 0.33025557, 0.62163083, 0.8369526 , + 0.52574254, 0.31510396, 0.06566058, 0.44301474, 0.27148304, + 0.30111694, 0.08374708, 0.11905168, 0.13817447, 0.35589242, + 0.20325078, 0.28584166, 0.13142756, 0.2040255 , 0.27938807, + 0.11800877, 0.03426762, 0.03599968, 0.01689927, 0.24260188, + -0.03815829, 0.19310196, 0.03499076, 0.10427147, 0.13123463, + 0.11359923, 0.2292006 , 0.51829568, 0.3945326 , 0.47134833, + 0.14234746, 0.17401509, 0.24885599, 0.29500217, 0.47867023, + 0.62890929, 0.34624069, 0.36416877, 0.3573325 , 0.69456097, + 1. , 0.45229606, 0.34654782, 0.12841021, 0.47623741, + 0.23457237, 0.06286334, 0.02768237, 0.00291156, 0.27546106, + 0.37390449, 0.29662638, 0.3605275 , 0.22705125, 0.26378275, + 0.22736693, 0.2490869 , 0.03850265, 0.0448428 , 0.07056435, + 0.13632041, 0.06124779, 0.13272672]) + +fc_48 = np.array([ 0.08348805, 0.25433363, 0.30718973, 0.21778039, -0.03452652, + -0.07191728, 0.04859293, 0.58070912, 0.3228884 , 0.03022416, + 0.096414 , 0.2067467 , 0.02934248, 0.11758467, 0.3395558 , + 0.2618392 , -0.04729804, 0.02403979, 0.33068409, 0.18542544, + 0.11157786, 0.24558548, 0.35921583, 0.14047045, 0.16072442, + 0.08725231, 0.04247384, 0.21984861, 0.15750137, 0.49704743, + 0.38107799, 0.18360744, 0.13072013, 0.03044336, 0.25038236, + 0.21535622, 0.04945266, 0.09180092, 0.07482232, 0.19016411, + -0.01817344, 0.30486544, 0.58977128, 0.33755094, 0.19382802, + 0.33391812, 0.38968696, 0.47185658, 1. , 0.74660754, + 0.25389967, 0.11339373, 0.30413138, 0.18101733, 0.22229945, + 0.48579973, 0.29879378, 0.14413303, 0.12710287, 0.61299693, + 0.3945326 , 0.03012108, 0.25349449, 0.36321952, 0.1842467 , + 0.09006558, 0.00976703, 0.01410019, 0.08918964, 0.24576006, + 0.72460302, 0.25165407, 0.31029377, 0.155622 , 0.07071528, + 0.26798855, 0.34361951, 0.15074631, 0.07711461, 0.13811093, + 0.20971265, 0.03274206, 0.08477825]) + +# Set the seed for the seed-based functional connectivity analysis. +seed = np.zeros_like(fc_48) +seed[60] = 1 +_cm = 'RdBu_r' + +# keywords +kws = {'edgecolors': 'k', + 'cmap': _cm, # 'Reds_r', #RdBu_r', # jet', + 'vmax': 0.5, + 'vmin': -.3, + 'alpha': None, + 'linewidth': 0.05} + +plot_surface_mpl_mv(vtx=vtx_l2k8_lhrh, + tri=tri_l2k8_lhrh, + hemi=hemi, + reorient='fs', + data=fc_60[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh], + shade_kwargs=kws) +"""_cm = 'BuGn' +kws = {'edgecolors': 'k', 'cmap': _cm, #$ 'Reds_r', #RdBu_r', # jet', + 'vmax': 1.0,'vmin': 0.0, 'alpha': 0.15, 'linewidth': 0.05} +mask = np.array([False]*40960) +mask[:vtx_l2k8_lhrh.shape[0]]= seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh] == 1 +mask1 = seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh] == 1 +plot_surface_mpl_mv(vtx=vtx_l2k8_lhrh,tri=tri_l2k8_lhrh,hemi=hemi,reorient='fs',data=seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh], + shade_kwargs=kws) """ + +fc_test_54 = np.array([4.04483037e-02, 7.32430791e-02, 7.76209493e-02, -4.46576574e-02, + 6.90320864e-02, 6.08216248e-02, 1.25866896e-01, 1.48259337e-01, + 1.46379082e-01, 1.60607085e-01, 1.99941445e-01, -4.73537989e-02, + 9.32955181e-02, 1.84146202e-01, 1.12089813e-01, 1.58828217e-01, + 1.37430776e-01, 3.38672182e-02, 8.48897343e-02, 1.80693377e-01, + 3.68792057e-02, -1.50098567e-02, 3.20046130e-02, 1.00762267e-02, + 1.56791288e-02, 2.76290084e-02, -3.96479495e-02, 2.80856785e-02, + 3.91929104e-02, 4.77790081e-02, 4.41821994e-02, 6.01223904e-02, + 5.10560110e-02, 7.76196128e-02, 1.12436866e-01, 2.73838335e-02, + -2.39388323e-02, 1.71982611e-02, 7.29186160e-02, 5.29034320e-02, + -1.99933514e-03, 7.06724821e-02, 5.93341413e-02, -4.13728002e-02, + 1.29869528e-04, 1.23697792e-01, 1.43169035e-01, 2.10767740e-01, + 3.11019690e-01, 1.38353247e-01, 1.86178990e-01, 2.70765146e-01, + -9.97468529e-03, 3.78946661e-01, 1.00000000e+00, 3.15421661e-01, + 1.52251462e-01, 1.33510163e-01, 1.12895918e-01, 1.14309818e-01, + 2.51496899e-01, 4.27383293e-02, 4.84373107e-02, 4.49691775e-02, + 6.54305540e-02, 6.93010133e-02, -3.49596325e-02, -1.16978796e-02, + -5.19326833e-02, 1.19297520e-01, 1.13479555e-01, 9.82750559e-02, + 1.37192609e-01, 1.10936182e-01, 1.36305410e-01, 1.35283865e-01, + 1.46130802e-01, 1.08109444e-01, 2.04031992e-02, -1.04834287e-02, + -5.35527088e-02, 3.55320361e-02, 5.05661223e-02]) + +fc_test_13 = np.array([ 0.02064446, -0.02125758, 0.06787312, -0.05500374, -0.03885027, + -0.01193319, 0.02590079, 0.07862877, 0.11140791, 0.12894951, + 0.29816439, -0.07257691, 0.39432721, 1. , 0.41285587, + 0.15905054, 0.02590495, 0.02414737, 0.16389443, 0.33678748, + 0.06684584, 0.02343256, 0.01415101, -0.00517585, 0.0300626 , + 0.14180168, 0.04087309, -0.02108286, 0.01249072, 0.02652301, + 0.0147603 , 0.05418977, 0.01977914, 0.03116899, 0.03394 , + 0.0515946 , -0.0622778 , -0.03677459, 0.02683279, 0.09717596, + 0.04846726, 0.0334019 , 0.06863676, -0.02671001, -0.04969303, + 0.04434678, 0.0240701 , 0.05507989, 0.09096978, 0.08247514, + 0.10845567, 0.22910861, -0.04563424, 0.13561862, 0.1841462 , + 0.16418864, 0.11336982, 0.11471379, 0.11565096, 0.05493549, + 0.27033029, 0.01242047, 0.03619441, 0.08456932, 0.05367578, + 0.10041786, -0.04090774, -0.01034624, -0.03077638, 0.08534603, + 0.06689338, 0.06317927, 0.07676879, 0.11862627, 0.05053353, + 0.10158778, 0.09397584, 0.05032458, 0.0745082 , 0.04716247, + 0.08466447, -0.05279958, 0.03736607]) + +fc_test_60 = np.array([3.51703550e-02, 5.18501633e-02, 3.00072240e-02, -3.22250993e-03, + 5.70771272e-03, 3.04967810e-02, 6.09221308e-02, 1.09301758e-01, + 1.10533317e-01, 1.28118958e-01, 2.10016752e-01, 1.09933921e-03, + 9.57229525e-02, 2.70330287e-01, 1.90774291e-01, 1.47761080e-01, + 8.97752877e-02, 1.11693599e-01, 1.95514360e-01, 2.72475220e-01, + 1.00020912e-01, 6.70112479e-02, 6.97584082e-02, 7.93095470e-02, + 8.93715930e-02, 7.42034261e-02, 3.39152899e-02, -4.49641276e-02, + -2.76792903e-02, 2.70008066e-02, -4.64044487e-03, 7.18338669e-02, + 2.04857838e-03, 8.14713997e-02, 1.17725733e-01, 7.42598533e-02, + 2.96991318e-03, -5.39451022e-02, 2.48211990e-02, 7.97761141e-02, + 1.85885648e-02, 1.56141833e-01, 3.34977809e-02, -3.00353457e-02, + -1.62362070e-02, 8.00848864e-02, 1.34135085e-01, 1.40980153e-01, + 1.77252780e-01, 1.22131104e-01, 2.10072865e-01, 2.10706170e-01, + -7.68060410e-04, 1.37570662e-01, 2.51496899e-01, 2.69339188e-01, + 2.11895067e-01, 1.73590178e-01, 2.42313753e-01, 2.07098915e-01, + 1.00000000e+00, 8.07081133e-02, 1.00150620e-01, 8.85262622e-02, + 9.18910872e-02, 1.83301361e-01, -1.60720267e-02, 3.92561061e-02, + 1.93154903e-03, 2.41607488e-01, 2.17235434e-01, 1.30866821e-01, + 2.47044743e-01, 2.56881848e-01, 1.50520001e-01, 2.11369632e-01, + 1.62973354e-01, 1.44820176e-01, 2.40272233e-02, 4.40490309e-02, + 3.92930927e-02, -4.07958819e-02, 4.43558957e-02]) + +fc_test_48 = np.array([ 0.04456317, 0.08430425, 0.10337379, -0.0751953 , 0.09872982, + 0.1291111 , 0.18314354, 0.20767963, 0.23303186, 0.19998463, + 0.14810478, -0.05109569, 0.05073159, 0.09096978, 0.08419494, + 0.17615717, 0.16614261, 0.0196335 , 0.09354732, 0.08995714, + 0.05359471, -0.00912246, 0.03514153, 0.02057001, 0.02612656, + 0.01342597, -0.0294798 , 0.00696833, 0.02931964, 0.03325134, + 0.07043438, 0.0502618 , 0.02825642, 0.18071165, 0.18340093, + 0.06189989, -0.00533773, 0.06097923, 0.10050126, 0.06412034, + 0.01178386, 0.17734297, 0.06325622, -0.02823774, -0.02826472, + 0.23305105, 0.33101713, 0.32282376, 1. , 0.31766683, + 0.28207603, 0.209288 , 0.00738726, 0.27140204, 0.31101969, + 0.22895553, 0.2234192 , 0.18360893, 0.16395947, 0.21453826, + 0.17725278, 0.0341468 , 0.03666224, 0.00832946, 0.05212808, + 0.07890769, 0.01180786, -0.01867186, -0.04819648, 0.11294647, + 0.16827467, 0.14320286, 0.16911658, 0.13737125, 0.19928486, + 0.21347781, 0.14709652, 0.16063771, 0.02247446, 0.04255669, + -0.00888587, 0.04339403, 0.07802365]) + +fc_test_41 = np.array([ 2.91574121e-02, 3.99662172e-02, 9.14379959e-02, -8.00334527e-03, + -2.50339264e-02, 3.63609413e-02, 1.06375335e-01, 1.06307706e-01, + 1.41449230e-01, 1.52185927e-01, 1.40082493e-01, 5.44236473e-03, + -2.84602123e-03, 3.34019030e-02, 7.24177217e-02, 1.29194284e-01, + 1.10487105e-01, 9.90362175e-02, 1.04988485e-01, 1.12635543e-01, + 1.02215036e-01, 3.63280972e-02, 9.13308055e-02, 6.25911165e-02, + 6.41443914e-02, 5.56257774e-02, 4.64617068e-02, -1.15959055e-02, + 9.60144121e-03, 5.03256381e-02, 3.93803611e-03, 5.98747396e-02, + -3.76220011e-02, 6.67289208e-02, 1.50632776e-01, 5.85915414e-02, + -5.79427852e-03, -2.67941385e-02, 2.60525511e-02, 8.43785417e-02, + 6.25090275e-02, 1.00000000e+00, 8.24955180e-02, -3.07733096e-02, + 3.37654868e-02, 1.69057126e-01, 1.98151129e-01, 2.99425215e-01, + 1.77342974e-01, 1.43791454e-01, 1.87007384e-01, 1.55162276e-01, + 2.97501823e-03, 1.24536132e-01, 7.06724821e-02, 1.09327160e-01, + 2.14722698e-01, 1.89543923e-01, 1.55433207e-01, 2.17070757e-01, + 1.56141833e-01, 7.86344152e-02, 8.39345856e-02, 9.96740939e-02, + 7.45054762e-02, 1.31522853e-01, 4.30021043e-02, -1.37114672e-02, + -3.84809861e-04, 1.87042770e-01, 2.42837501e-01, 1.73869281e-01, + 2.42828040e-01, 1.65126795e-01, 3.06438319e-01, 2.48172135e-01, + 2.74101333e-01, 2.87125048e-01, 3.19166679e-03, 1.65886750e-02, + 3.04314795e-02, 3.21617761e-02, 4.86736314e-03]) + +seed = np.zeros_like(fc_48) +seed[48] = 1 +_cm = 'RdBu_r' +kws = {'edgecolors': 'k', + 'cmap': _cm, # 'Reds_r', #RdBu_r', # jet', + 'vmax': 0.35, + 'vmin': -.2, + 'alpha': None, + 'linewidth': 0.05} + +plot_surface_mpl_mv(vtx=vtx_l2k8_lhrh, + tri=tri_l2k8_lhrh, + hemi=hemi, + reorient='fs', + data=fc_test_48[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh], + shade_kwargs=kws) +"""_cm = 'BuGn' +kws = {'edgecolors': 'k', 'cmap': _cm, #$ 'Reds_r', #RdBu_r', # jet', + 'vmax': 1.0,'vmin': 0.0, 'alpha': 0.15, 'linewidth': 0.05} +mask = np.array([False]*40960) +mask[:vtx_l2k8_lhrh.shape[0]]= seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh] == 1 +mask1 = seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh] == 1 +plot_surface_mpl_mv(vtx=vtx_l2k8_lhrh,tri=tri_l2k8_lhrh,hemi=hemi,reorient='fs',data=seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh], + shade_kwargs=kws)""" + +end_time = time.time() +print('Running time is {0} \'s'.format(end_time - start_time )) + +