# 导入工具包
import h5py
import rcca
import sys
import numpy as np
import cortex
zscore = lambda d: (d-d.mean(0))/d.std(0)

data = []
vdata = []
numSubjects = 3
# subjects 是3个受试者列表.
subjects = ['S1', 'S2', 'S3']
# xfms 是Pycortex中变换名称的列表，该名称用于对齐每个受试者的功能和解剖数据。
xfms = ['S1_xfm', 'S2_xfm', 'S3_xfm']
dataPath ="./data/VoxelResponses_subject%d.mat"
for subj in range(numSubjects):# 打开数据f = h5py.File(dataPath % (subj+1),'r')#  获取数据大小datasize = (int(f["ei"]["datasize"].value[2]),int(f["ei"]["datasize"].value[1]),int(f["ei"]["datasize"].value[0]))# 从Pycortex获取皮质面罩mask = cortex.db.get_mask(subjects[subj], xfms[subj], type = 'thick')# 获取该受试者的训练数据data_subj = np.nan_to_num(zscore(np.nan_to_num(f["rt"].value.T)))data.append(data_subj.reshape((data_subj.shape[0],)+datasize)[:, mask])# 获取受试者的验证数据vdata_subj = np.nan_to_num(zscore(np.nan_to_num(f["rv"].value.T)))vdata.append(vdata_subj.reshape((vdata_subj.shape[0],)+datasize)[:, mask])

# 这里设置1e-4和1e2之间的一系列正则化值
regs = np.array(np.logspace(-4, 2, 10))# 这里考虑3到10之间的成分数量
numCCs = np.arange(3, 11)# 初始化cca模型
cca = rcca.CCACrossValidate(numCCs=numCCs, regs=regs)

"""
说明：
由于数据量大，此分析的计算量很大。
在笔记本中运行它会花费大量时间，因此建议对其进行并行化和/或在计算机群集上运行它，然后加载结果以进行可视化。
"""# 利用cca训练数据
cca.train(data)# 利用cca对验证数据进行验证
cca.validate(vdata)# 计算方差，解释每个体素中的验证响应
cca.compute_ev(vdata)# 保存分析结果
cca.save("./data/CCA_results.hdf5")

# 导入可视化工具包
%matplotlib inline
import matplotlib.pyplot as plt# 导入Brewer色彩图以进行可视化
from brewer2mpl import qualitativenSubj = len(cca.corrs)
nBins = 30
bmap = qualitative.Set1[nSubj]
f = plt.figure(figsize = (8, 6))
ax = f.add_subplot(111)
for s in range(nSubj):# 绘制所有三个对象的所有体素之间的相关性直方图ax.hist(cca.corrs[s], bins = nBins, color = bmap.mpl_colors[s], histtype="stepfilled", alpha = 0.6)
plt.legend(['Subject 1', 'Subject 2', 'Subject 3'], fontsize = 16)
ax.set_xlabel('Prediction correlation', fontsize = 20)
ax.set_ylabel('Number of voxels', fontsize = 20)
ax.set_title("Prediction performance across voxels", fontsize = 20)
# p <0.05时的显着性阈值(针对多次比较进行了校正)
# 重要性是使用渐近方法计算的（有关详细信息，请参见论文文本）
thresh = 0.0893
ax.axvline(x = thresh, ymin = 0, ymax = 7000, linewidth = 2, color = 'k')
ax.text(thresh+0.05, 5000, 'p<0.05', fontsize = 16)
ax.set_xticklabels(0.1*np.arange(-8, 11, 2), fontsize = 16)
ax.set_yticklabels(np.arange(0, 10000, 1000), fontsize = 16)

import cortex
from matplotlib import cm
from copy import deepcopy
subj = 0
subjName = "S1"
subjTransform = "S1_xfm"
corrs = deepcopy(cca.corrs[subj])
# 将所有低于显着性阈值的体素设置为0
corrs[corrs<thresh] = 0
_ = cortex.quickflat.make_figure(cortex.Volume(corrs, subjName, subjTransform, cmap = cm.PuBuGn_r,vmin = 0., vmax = 1.), with_curvature = True)