APR6p_patterns_plot_source_stats.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8640b1bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import mne\n",
    "#%gui qt\n",
    "#import matplotlib\n",
    "\n",
    "#%matplotlib qt\n",
    "from nilearn.plotting import plot_stat_map\n",
    "from nilearn.image import index_img\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from stormdb.access import Query\n",
    "from pickle import load\n",
    "from scipy import stats as spstats\n",
    "from mne.datasets import sample\n",
    "from mne.stats import spatio_temporal_cluster_1samp_test\n",
    "import os\n",
    "import pickle\n",
    "import pandas as pd\n",
    "from copy import deepcopy\n",
    "import warnings\n",
    "import os.path as op\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning) \n",
    "# os.environ['ETS_TOOLKIT'] = 'qt4'\n",
    "# os.environ['QT_API'] = 'pyqt5'\n",
    "# %gui qt\n",
    "#mne.viz.set_3d_backend(\"notebook\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7509def",
   "metadata": {},
   "outputs": [],
   "source": [
    "proj_name = 'MINDLAB2020_MEG-AuditoryPatternRecognition'\n",
    "wdir = '/projects/' + proj_name + '/scratch/working_memory/'\n",
    "stats_dir = wdir + 'results/stats/'\n",
    "data_dir = wdir + 'averages/data/'\n",
    "figures_dir = wdir + 'results/figures/'\n",
    "\n",
    "subs_dir = '/projects/' + proj_name + '/scratch/fs_subjects_dir/'\n",
    "sample_path = sample.data_path()\n",
    "sample_subjects_dir = sample_path + '/subjects'\n",
    "src_sample = mne.read_source_spaces(subs_dir +\n",
    "                                '/fsaverage/bem/fsaverage-vol-5-src.fif')\n",
    "# src_sample = mne.read_source_spaces(subs_dir +\n",
    "#                                 'fsaverage/bem/fsaverage_ico4_vol-src.fif')\n",
    "#label_file = subs_dir + 'fsaverage/mri/aparc.DKTatlas+aseg.mgz'\n",
    "#label_file = sample_subjects_dir + '/fsaverage/mri/aparc.a2009s+aseg.mgz'\n",
    "label_file = sample_subjects_dir + '/fsaverage/mri/aparc+aseg.mgz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7801816c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77fc7c09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read stats\n",
    "periods = {#'all':[0, 4],\n",
    "           #'list_vs_im': [0,4],\n",
    "           'L1': [.2,.5],\n",
    "           'L2': [.7,1.],\n",
    "           'L3': [1.2,1.5],\n",
    "           'I1': [2,4],\n",
    "           'I2': [2.5,3.5]\n",
    "           }#,'retrieval': [4,6.25]}{'all': [0.,4]}#\n",
    "stats = {}\n",
    "\n",
    "for p in periods:\n",
    "    csfname = '{}patterns_source_stats_{}_peaks.p'.format(stats_dir, p)\n",
    "    #csfname = '{}patterns_source_stats_{}.p'.format(stats_dir, p)\n",
    "    cfile = open(csfname,'rb')\n",
    "    stats[p] = pickle.load(cfile)\n",
    "    cfile.close()\n",
    "print(stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fd29e1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "for c in stats:\n",
    "    print(c)\n",
    "    for s in stats[c]:\n",
    "        print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3fb314c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load and morph a source time course\n",
    "dfname = data_dir + '0021_LZW/0021_LZW_patterns_sources_task_sensor_lf_0.05_hf_None_tstep_0.025_twin_0.05_localized.p'\n",
    "dfile = open(dfname,'rb')\n",
    "all_stc = pickle.load(dfile)\n",
    "dfile.close()\n",
    "stc = deepcopy(all_stc['maintenance1'])\n",
    "del all_stc\n",
    "\n",
    "# load source morph\n",
    "morph = mne.read_source_morph(subs_dir + '0021_LZW/bem/0021_LZW_vol-morph.h5')\n",
    "stc = morph.apply(stc)\n",
    "print(stc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "790aa661",
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = mne.get_volume_labels_from_aseg(label_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "925b2591",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c86f243e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#labels = mne.read_labels_from_annot('fsaverage','aparc.a2009s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aa5030b",
   "metadata": {},
   "outputs": [],
   "source": [
    "labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15ebf880",
   "metadata": {},
   "outputs": [],
   "source": [
    "clabels = ['Left-Thalamus-Proper',\n",
    "           'Right-Thalamus-Proper',\n",
    "           'Left-VentralDC',\n",
    "           'Right-VentralDC',\n",
    "           'Right-Hippocampus',\n",
    "           'Left-Hippocampus',\n",
    "           'ctx_lh_G_oc-temp_med-Parahip',\n",
    "           'ctx_rh_G_oc-temp_med-Parahip',\n",
    "           'ctx_lh_G_and_S_cingul-Ant',\n",
    "           'ctx_rh_G_and_S_cingul-Ant',\n",
    "           'ctx_lh_G_and_S_cingul-Mid-Ant',\n",
    "           'ctx_rh_G_and_S_cingul-Mid-Ant',\n",
    "           'ctx_lh_G_and_S_cingul-Mid-Post',\n",
    "           'ctx_rh_G_and_S_cingul-Mid-Post',\n",
    "           'ctx_lh_G_temp_sup-Plan_tempo',\n",
    "           'ctx_rh_G_temp_sup-Plan_tempo',\n",
    "           'ctx_lh_G_temp_sup-Lateral',\n",
    "           'ctx_rh_G_temp_sup-Lateral',\n",
    "           'Left-Pallidum',\n",
    "           'Right-Pallidum',\n",
    "           'Left-Caudate',\n",
    "           'Right-Caudate',\n",
    "           'ctx_lh_S_orbital-H_Shaped',\n",
    "           'ctx_rh_S_orbital-H_Shaped',\n",
    "           'ctx_lh_G_front_inf-Opercular',\n",
    "           'ctx_rh_G_front_inf-Opercular',\n",
    "           'ctx_lh_G_front_inf-Orbital',\n",
    "           'ctx_rh_G_front_inf-Orbital',\n",
    "           'ctx_lh_G_front_inf-Triangul',\n",
    "           'ctx_rh_G_front_inf-Triangul',\n",
    "           'ctx_lh_G_front_middle',\n",
    "           'ctx_rh_G_front_middle',\n",
    "           'ctx_lh_G_front_sup',\n",
    "           'ctx_rh_G_front_sup',\n",
    "           'ctx_lh_Lat_Fis-ant-Horizont',\n",
    "           'ctx_rh_Lat_Fis-ant-Horizont']\n",
    "\n",
    "stc_labels = {}\n",
    "se_labels = {}\n",
    "stc_mask = {}\n",
    "conds = ['maintenance1','manipulation1','interaction']\n",
    "for c in conds:\n",
    "    cdata = stc.copy().crop(-.05,4)\n",
    "    cmask = stc.copy().crop(-.05,4)\n",
    "    cse = stc.copy().crop(-.05,4)\n",
    "    emask = stats['encoding'][c]['mask'].copy()\n",
    "    dmask = stats['delay'][c]['mask'].copy()\n",
    "\n",
    "    cdata.data = np.hstack((stats['encoding'][c]['data_mean'],stats['delay'][c]['data_mean'])) # * stats_results[c]['mask']\n",
    "    cmask.data = np.hstack((emask,dmask))\n",
    "    print(cmask.data)\n",
    "    cse.data = np.hstack((stats['encoding'][c]['data_sd'],stats['delay'][c]['data_sd'])) / np.sqrt(stats['encoding'][c]['n'])# * stats_results[c]['mask'] \n",
    "#    cdata.data = stats['encoding'][c]['data_mean'].T # * stats_results[c]['mask'] \n",
    "#    cse.data = stats['encoding'][c]['data_sd'].T / np.sqrt(stats['encoding'][c]['n'])# * stats_results[c]['mask'] \n",
    "    stc_labels[c] = cdata.extract_label_time_course(labels = [label_file,clabels], src = src_sample, mode = 'auto')\n",
    "    stc_mask[c] = cmask.extract_label_time_course(labels = [label_file,clabels], src = src_sample, mode = 'auto')\n",
    "    se_labels[c] = cse.extract_label_time_course(labels = [label_file,clabels], src = src_sample, mode = 'auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a6b7a6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "times = stc.copy().crop(-.025,4).times \n",
    "periods = {'encoding': [0,2], 'delay': [2.0,3.975]}\n",
    "blocks = ['maintenance1','manipulation1']\n",
    "cur_labs = ['Right Planum Temporale','Right Hippocampus','Right Pallidum',\n",
    "            'Left Thalamus','Right Ventral Diencephalon',\n",
    "            'Right Caudate','Right Orbital Sulcus','Anterior Lateral Fissure',\n",
    "            'Posterior Cingulate']\n",
    "cur_labs_codes = ['ctx_rh_G_temp_sup-Plan_tempo','Right-Hippocampus','Right-Pallidum',\n",
    "                  'Left-Thalamus-Proper','Right-VentralDC',\n",
    "                  'Right-Caudate','ctx_rh_S_orbital-H_Shaped','ctx_lh_Lat_Fis-ant-Horizont',\n",
    "                  'ctx_rh_G_and_S_cingul-Mid-Post']\n",
    "lab_ix = []\n",
    "for ll in cur_labs_codes:\n",
    "    lab_ix += [np.where(np.array(clabels) == ll)[0][0]]\n",
    "    print(lab_ix)\n",
    "#lab_ix = [11,15]\n",
    "titles = ['maintenance','manipulation']\n",
    "sidx = -1\n",
    "fig, axes = plt.subplots(ncols=2,nrows=len(cur_labs), figsize = (12.5,len(cur_labs)*5))\n",
    "\n",
    "for clix, cl in enumerate(cur_labs):\n",
    "    for bidx,b in enumerate(blocks):\n",
    "        lidx = lab_ix[clix]\n",
    "        sidx = sidx + 1\n",
    "        rix, cix = sidx//2,sidx%2\n",
    "        #axes[rix, cix].plot(times,stc_labels[b + '/' + 'mel1'][lidx],color='b')\n",
    "        #axes[rix, cix].plot(times,stc_labels[b + '/' + 'mel2'][lidx],color='r')\n",
    "        ci_upper = stc_labels[b][lidx] + 1.96*se_labels[b][lidx]\n",
    "        ci_lower = stc_labels[b][lidx] - 1.96*se_labels[b][lidx]\n",
    "        ccmask = np.array(stc_mask[b][lidx] != 0).astype(float)\n",
    "        ccmask[ccmask==0.] = np.nan\n",
    "        #axes[rix, cix].fill_between(times, ci_lower, ci_upper, color='k', alpha=.2)\n",
    "        axes[rix, cix].plot(times,stc_labels[b][lidx],color='k')\n",
    "        axes[rix, cix].plot(times,stc_labels[b][lidx]*ccmask,color='k',linewidth=4)\n",
    "        axes[rix, cix].set_title(titles[bidx] + ' ' + cl)\n",
    "        axes[rix, cix].set_ylim(-3,3)\n",
    "        axes[rix, cix].set_xlim(-.1,4)\n",
    "        axes[rix, cix].set_xlabel('time (s)')\n",
    "        axes[rix, cix].set_ylabel('source activation (a.u.)')\n",
    "        axes[rix, cix].axhline(0., color='k')\n",
    "        axes[rix, cix].axvline(0., color='k')\n",
    "        axes[rix, cix].axvline(2., color='k',linestyle='--')\n",
    "        axes[rix, cix].axvline(.5, color='k',linestyle=':')\n",
    "        axes[rix, cix].axvline(1, color='k',linestyle=':')\n",
    "        #axes[rix, cix].legend([clabels[cls] for cls in [0,8]])\n",
    "        #axes[rix, cix].legend(['patterns'])#,'difference'])\n",
    "plt.tight_layout()\n",
    "plt.savefig(figures_dir + 'patterns_mels_sources_labels_diff.pdf',orientation='portrait')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc7c3be6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_cluster_peaks(cstat,p_threshold=.05):\n",
    "        sigix = np.where(cstat['pvals']<=p_threshold)[0]\n",
    "        vx, vy, masks, p = [], [], [], []\n",
    "        for six in sigix:\n",
    "            cmask = cstat['clusters'][six].copy().astype(int).T\n",
    "            cdata = cstat['tvals']*cmask\n",
    "            cvx, cvy = np.unravel_index(np.abs(cdata).argmax(), cdata.shape)\n",
    "            vx += [cvx]\n",
    "            vy += [cvy]\n",
    "            masks += [cmask]\n",
    "            p += [cstat['pvals'][six]]\n",
    "        vx = np.array(vx)\n",
    "        vy = np.array(vy)\n",
    "        return vx, vy, masks, p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06cc479b",
   "metadata": {},
   "outputs": [],
   "source": [
    "for pix, pval in enumerate(stats[p][c]['pvals']):\n",
    "    print(pval)\n",
    "    print(np.min(np.abs(stats[p][c]['tvals']))*stats[p][c]['clusters'][pix].astype(int))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "548ae865",
   "metadata": {},
   "outputs": [],
   "source": [
    "conds = ['maintenance1','manipulation1','interaction']#'maintenance1','manipulation1',]\n",
    "cnames = ['recall','manipulation','difference']\n",
    "fig,ax = plt.subplots(nrows=1,ncols=1,figsize=(10,10))\n",
    "for p in stats:\n",
    "    for cix,c in enumerate(conds):#stats[p]:\n",
    "        vx, vy, masks, _ = get_cluster_peaks(stats[p][c])\n",
    "        for v in range(vx.shape[0]):\n",
    "            cdata = stats[p][c]['data_mean']*masks[v]\n",
    "            cond_name = cnames[cix]\n",
    "            cluster_name = p + ' (' + cond_name + ' - cluster ' + str(v+1) + ')'\n",
    "            ch = ax.plot(stc.copy().crop(tmin=periods[p][0],tmax=periods[p][1]).times,\n",
    "                     cdata.mean(axis=0),label=cluster_name)\n",
    "            plt.ylim(-.2,.2)\n",
    "            plt.xlim(-.1,4)\n",
    "            plt.xlabel('time (s)')\n",
    "            plt.ylabel('source activation (a.u.)')\n",
    "            plt.axhline(0., color='k')\n",
    "            plt.axvline(0., color='k')\n",
    "            plt.axvline(2., color='k',linestyle='--')\n",
    "            plt.axvline(.5, color='k',linestyle=':')\n",
    "            plt.axvline(1, color='k',linestyle=':')\n",
    "            plt.legend()\n",
    "plt.tight_layout()\n",
    "#plt.savefig(figures_dir + 'patterns_mels_sources_cluster.pdf',orientation='portrait')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebb1ede5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#periods = {'encoding': [0,2], 'delay': [2.0,3.975]}\n",
    "#periods = {'L2': [1.3,1]}#, 'delay': [2.0,3.975]}\n",
    "cperiods = ['L1','L2','L3','I1']#,'I2']#'list_vs_im'\n",
    "conds = ['maintenance1','manipulation1','interaction']\n",
    "label_sel = [l for l in labels if ('White' not in l) & ('choroid' not in l) & ('vessel' not in l) & ('callosum' not in l) &\n",
    "              ('Vent' not in l) & ('wall' not in l) & ('Unknown' not in l) & ('unknown' not in l) & ('Cerebellum' not in l)\n",
    "              & ('CC_' not in l) & ('CSF' not in l) & ('Optic-Chiasm' not in l) & ('WM-hypointensities' not in l)]\n",
    "MA = {}\n",
    "clust_report = {'period': [], 'condition': [], 'cluster No': [],\n",
    "                 'x': [], 'y': [], 'z': [],'p-value': [],\n",
    "                 'peak t-value': [], 'peak activation': []} \n",
    "for p in cperiods:\n",
    "    print(p)\n",
    "    MA[p] = {}\n",
    "    conds = stats[p].keys()\n",
    "    for c in conds:#stats[p]: \n",
    "        print(c)\n",
    "        cstc = stc.copy().crop(tmin=periods[p][0],tmax=periods[p][1])\n",
    "        vx, vy, masks, pvals = get_cluster_peaks(stats[p][c])\n",
    "        MA[p][c] = np.zeros((vx.shape[0],len(label_sel)))\n",
    "        for v in range(vx.shape[0]):\n",
    "            ctime = cstc.times[vy[v]]\n",
    "            print('cluster',v+1,' - peak: ', np.round(ctime,3), 's')\n",
    "            print(np.unique(masks[v]))\n",
    "            cstc.data = stats[p][c]['tvals'].copy()*masks[v]\n",
    "            csign = np.sign(cstc.data.mean())\n",
    "            cstc.data = np.abs(cstc.data)\n",
    "            pos, _ = cstc.get_peak()\n",
    "            pix, _ = cstc.get_peak(vert_as_index=True)\n",
    "            coords = np.round(src_sample[0]['rr'][pos]*1000)\n",
    "            clust_report['period'] += [p]\n",
    "            clust_report['condition'] += [c]\n",
    "            clust_report['cluster No'] += [v+1]\n",
    "            clust_report['x'] += [coords[0]]\n",
    "            clust_report['y'] += [coords[1]]\n",
    "            clust_report['z'] += [coords[2]]\n",
    "            clust_report['p-value'] += [pvals[v]]\n",
    "            clust_report['peak t-value'] += [np.round(cstc.data[pix,0]*csign,2)]\n",
    "            clust_report['peak activation'] += [np.round(stats[p][c]['data_mean'][pix,0],2)]\n",
    "\n",
    "            # MA[p][c][v] = np.squeeze(mne.extract_label_time_course(cstc,[label_file,label_sel],\n",
    "            #                                                         src=src_sample,mode='max'))*csign\n",
    "clust_report = pd.DataFrame(clust_report)         \n",
    "print(clust_report)\n",
    "clust_report.to_csv(op.join(stats_dir,'patterns_source_stats_report.csv'),index=False)\n",
    "\n",
    "            #MA[p][c][v] = np.max(np.abs(mne.extract_label_time_course(cstc,[label_file,labels],src=src_sample)),axis=1)*np.sign(cstc.data.mean())#.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e3c9874",
   "metadata": {},
   "outputs": [],
   "source": [
    "hems = {'right': ['rh','Right'],'left': ['lh','Left']}\n",
    "MA_df = {'cluster': [], 'label': [], 'hemisphere': [], 'tval': []}\n",
    "for p in MA:\n",
    "    for c in MA[p]:#stats[p]:\n",
    "        for cl in range(MA[p][c].shape[0]):\n",
    "            act_ix = np.argsort(-np.abs(MA[p][c][cl,:]))\n",
    "            #print( '\\n', p, c,' cluster ', cl + 1)\n",
    "            for h in hems:\n",
    "                #print(h)\n",
    "                for ix in act_ix:\n",
    "                    if (np.abs(MA[p][c][cl,ix]) >= 2) & (hems[h][0] in label_sel[ix] or hems[h][1] in label_sel[ix]): \n",
    "                        print(label_sel[ix],np.round(MA[p][c][cl,ix],2))\n",
    "                        MA_df['cluster'] += ['{} {} cluster {}'.format(p, c, cl + 1)]\n",
    "                        MA_df['label'] += [label_sel[ix]]\n",
    "                        MA_df['hemisphere'] += [h]\n",
    "                        MA_df['tval'] += [np.round(MA[p][c][cl,ix],2)]\n",
    "                             \n",
    "MA_df = pd.DataFrame(MA_df)\n",
    "print(MA_df)     \n",
    "MA_df.to_csv(op.join(stats_dir,'patterns_stats_report_aparc.csv'),index=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "376bf7a2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3acc33c",
   "metadata": {},
   "outputs": [],
   "source": [
    "conds = ['maintenance1','manipulation1','interaction']#'maintenance1','manipulation1',]\n",
    "for p in stats:\n",
    "    for c in conds:#stats[p]:\n",
    "        \n",
    "        cstc = stc.copy().crop(tmin=periods[p][0],tmax=periods[p][1])\n",
    "        vx, vy, masks, _ = get_cluster_peaks(stats[p][c])\n",
    "\n",
    "        for v in range(vy.shape[0]):\n",
    "            cstc.data = stats[p][c]['data_mean']*masks[v]\n",
    "            ctime = cstc.times[vy[v]]\n",
    "            #cstc = cstc.copy().crop(tmin=ctime-.1,tmax=ctime+.1).mean()\n",
    "            cstc = cstc.copy().mean()\n",
    "\n",
    "            img = cstc.as_volume(src_sample,  mri_resolution=False)\n",
    "            cname = '{} {} cluster {} - {}s'.format(p,c,v+1,np.round(ctime,2))\n",
    "            print(cname)\n",
    "            cfig = plt.figure()\n",
    "            brain = plot_stat_map(img, subs_dir + 'fsaverage/mri/T1.mgz',draw_cross=False,black_bg=True,\n",
    "                          #output_file = figures_dir + 'patterns_sources_' + cname + '.pdf',\n",
    "                          vmax=1,symmetric_cbar=True,figure=cfig, threshold=.025)\n",
    "            plt.show()\n",
    "            \n",
    "            # brain = cstc.copy().crop(tmin=ctime-.1,tmax=ctime+.1).mean().plot(subjects_dir=subs_dir,#size=1600,\n",
    "            #                                                              clim={'kind': 'value', 'pos_lims': [0.001,.2,1.2]},#time_viewer=True,\n",
    "            #                      subject='fsaverage', mode = 'stat_map',src=src_sample, initial_time = ctime)#,)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0d6ab9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import datasets\n",
    "from nilearn import surface\n",
    "from nilearn import plotting\n",
    "fsaverage = datasets.fetch_surf_fsaverage()\n",
    "fsaverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c2e4a62",
   "metadata": {},
   "outputs": [],
   "source": [
    "#conds = ['maintenance1','manipulation1','interaction']#'maintenance1','manipulation1',]\n",
    "#csfname = stats[].get_keys()\n",
    "views = ['lateral','medial']\n",
    "hemi = ['left','right']\n",
    "surfs = [fsaverage.pial_left, fsaverage.pial_right]\n",
    "surfs2 = [fsaverage.infl_left, fsaverage.infl_right]\n",
    "bgs = [fsaverage.sulc_left, fsaverage.sulc_right]\n",
    "for p in stats:\n",
    "    conds = stats[p].keys()\n",
    "    for c in conds:#stats[p]:   \n",
    "        vx, vy, masks, pvals = get_cluster_peaks(stats[p][c],p_threshold=.1)\n",
    "        print(vy)\n",
    "        for v in range(vy.shape[0]):\n",
    "            cstc = stc.copy().crop(tmin=periods[p][0],tmax=periods[p][1])\n",
    "            cstc.data = stats[p][c]['data_mean']*masks[v]\n",
    "            ctime = cstc.times[vy[v]].copy()\n",
    "\n",
    "            #cstc = cstc.copy().crop(tmin=ctime-.1,tmax=ctime+.1).mean()\n",
    "            cstc = cstc.copy().mean()\n",
    "            img = cstc.as_volume(src_sample,  mri_resolution=False)\n",
    "            \n",
    "            cname = '{} {} cluster {} - {}s (p = {})'.format(p,c,v+1,np.round(ctime,2),np.round(pvals[v],3))\n",
    "            print(cname)\n",
    "            cfig, caxis = plt.subplots(nrows=2,ncols=2,subplot_kw={'projection': '3d'})\n",
    "            for vix, vi in enumerate(views):\n",
    "                for hix,h in enumerate(hemi):\n",
    "                    surf = surface.vol_to_surf(img, surfs[hix])\n",
    "                    plotting.plot_surf_stat_map(surfs2[hix],surf, hemi=h,axes=caxis[vix,hix],\n",
    "                                                   symmetric_cbar=True,view=vi, vmax=1,threshold=.025,\n",
    "                                                    colorbar=True, bg_map=bgs[hix],darkness=1)#,#,\n",
    "                                                   # engine='matplotlib',bg_on_data=True)\n",
    "            plt.show()\n",
    "            #cfig.savefig(figures_dir + 'patterns_sources_' + cname + '_surface.pdf')\n",
    "            # brain = cstc.copy().crop(tmin=ctime-.1,tmax=ctime+.1).mean().plot(subjects_dir=subs_dir,#size=1600,\n",
    "            #                                                              clim={'kind': 'value', 'pos_lims': [0.001,.2,1.2]},#time_viewer=True,\n",
    "            #                      subject='fsaverage', mode = 'stat_map',src=src_sample, initial_time = ctime)#,)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ea05fea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#conds = ['maintenance1','manipulation1','interaction']#'maintenance1','manipulation1',]\n",
    "views = ['lateral','medial']\n",
    "hemi = ['left','right']\n",
    "surfs = [fsaverage.pial_left, fsaverage.pial_right]\n",
    "surfs2 = [fsaverage.infl_left, fsaverage.infl_right]\n",
    "bgs = [fsaverage.sulc_left, fsaverage.sulc_right]\n",
    "for p in stats:\n",
    "    conds = stats[p].keys()\n",
    "    for c in conds:#stats[p]:   \n",
    "        # vx, vy, masks, pvals = get_cluster_peaks(stats[p][c])\n",
    "        # print(vy)\n",
    "    #for v in range(vy.shape[0]):\n",
    "        cstc = stc.copy().crop(tmin=periods[p][0],tmax=periods[p][1])\n",
    "        cstc.data = stats[p][c]['data_mean']*stats[p][c]['mask']\n",
    "        #ctime = cstc.times[vy[v]].copy()\n",
    "\n",
    "        #cstc = cstc.copy().crop(tmin=ctime-.1,tmax=ctime+.1).mean()\n",
    "        cstc = cstc.copy().mean()\n",
    "        img = cstc.as_volume(src_sample,  mri_resolution=False)\n",
    "        \n",
    "        cname = '{} {}'.format(p,c)\n",
    "        print(cname)\n",
    "        cfig, caxis = plt.subplots(nrows=2,ncols=2,subplot_kw={'projection': '3d'})\n",
    "        for vix, vi in enumerate(views):\n",
    "            for hix,h in enumerate(hemi):\n",
    "                surf = surface.vol_to_surf(img, surfs[hix])\n",
    "                plotting.plot_surf_stat_map(surfs2[hix],surf, hemi=h,axes=caxis[vix,hix],\n",
    "                                                symmetric_cbar=True,view=vi, vmax=1,threshold=.025,\n",
    "                                                colorbar=True, bg_map=bgs[hix],darkness=1)#,#,\n",
    "                caxis[vix,hix].set_rasterized(True) #To save more efficiently\n",
    "\n",
    "                                                # engine='matplotlib',bg_on_data=True)\n",
    "        plt.show()\n",
    "        cfig.savefig(figures_dir + 'patterns_sources_' + cname + '_all_clust_surface.pdf',dpi=300)\n",
    "            # brain = cstc.copy().crop(tmin=ctime-.1,tmax=ctime+.1).mean().plot(subjects_dir=subs_dir,#size=1600,\n",
    "            #                                                              clim={'kind': 'value', 'pos_lims': [0.001,.2,1.2]},#time_viewer=True,\n",
    "            #                      subject='fsaverage', mode = 'stat_map',src=src_sample, initial_time = ctime)#,)\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c763f618",
   "metadata": {},
   "outputs": [],
   "source": [
    "conds = ['interaction']#'maintenance1','manipulation1',]\n",
    "views = ['lateral','medial']\n",
    "hemi = ['left','right']\n",
    "surfs = [fsaverage.pial_left, fsaverage.pial_right]\n",
    "surfs2 = [fsaverage.infl_left, fsaverage.infl_right]\n",
    "bgs = [fsaverage.sulc_left, fsaverage.sulc_right]\n",
    "for p in ['I2']:\n",
    "    for c in conds:#stats[p]:   \n",
    "        # vx, vy, masks, pvals = get_cluster_peaks(stats[p][c])\n",
    "        # print(vy)\n",
    "    #for v in range(vy.shape[0]):\n",
    "        cstc = stc.copy().crop(tmin=periods[p][0],tmax=periods[p][1])\n",
    "        cstc.data = stats[p][c]['data_mean']*(np.abs(stats[p][c]['tvals']) > 2)\n",
    "        #ctime = cstc.times[vy[v]].copy()\n",
    "\n",
    "        #cstc = cstc.copy().crop(tmin=ctime-.1,tmax=ctime+.1).mean()\n",
    "        cstc = cstc.copy().mean()\n",
    "        img = cstc.as_volume(src_sample,  mri_resolution=False)\n",
    "        \n",
    "        cname = '{} {}'.format(p,c)\n",
    "        print(cname)\n",
    "        cfig, caxis = plt.subplots(nrows=2,ncols=2,subplot_kw={'projection': '3d'})\n",
    "        for vix, vi in enumerate(views):\n",
    "            for hix,h in enumerate(hemi):\n",
    "                surf = surface.vol_to_surf(img, surfs[hix])\n",
    "                plotting.plot_surf_stat_map(surfs2[hix],surf, hemi=h,axes=caxis[vix,hix],\n",
    "                                                symmetric_cbar=True,view=vi, vmax=1,threshold=.025,\n",
    "                                                colorbar=True, bg_map=bgs[hix],darkness=1)#,#,\n",
    "                                                # engine='matplotlib',bg_on_data=True)\n",
    "                caxis[vix,hix].set_rasterized(True) #To save more efficiently\n",
    "\n",
    "                                                # engine='matplotlib',bg_on_data=True)\n",
    "        plt.show()\n",
    "        cfig.savefig(figures_dir + 'patterns_sources_' + cname + '_thresholded_surface.pdf',dpi=300)\n",
    "        #cfig.savefig(figures_dir + 'patterns_sources_' + cname + '_all_clust_surface.pdf')\n",
    "            # brain = cstc.copy().crop(tmin=ctime-.1,tmax=ctime+.1).mean().plot(subjects_dir=subs_dir,#size=1600,\n",
    "            #                                                              clim={'kind': 'value', 'pos_lims': [0.001,.2,1.2]},#time_viewer=True,\n",
    "            #                      subject='fsaverage', mode = 'stat_map',src=src_sample, initial_time = ctime)#,)\n",
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}