{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project pRF analysis results to pycortex surfaces\n",
    "Make sure you have created the pycortex entries for these individuals (https://github.com/VisionandCognition/NHP-pycortex)        \n",
    "NB! the pycortex database for these subjects is available in `<..>/Preprocessed_data/MRI/pycortex_db`       \n",
    "\n",
    "Strategy:    \n",
    "- Load all (unthresholded) results into sorted numpy arrays with nibabel\n",
    "- Also create pycortex volume objects\n",
    "- Save them in an ordered dictionary\n",
    "- Define some functions for masking\n",
    "- Perform any additional post-processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import the required modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cortex\n",
    "import nibabel as nib \n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os, shutil, copy\n",
    "\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set subject name and path to FitResults    \n",
    "NB! adjust the paths to where you saved things"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subs = ['M01', 'M02']\n",
    "\n",
    "FitResPath = os.path.join(<add path to this folder>,'Analysis_code','FitResults')\n",
    "ManualMaskPath = os.path.join(<add path to this folder>,'Analysis_code','Preprocessed_data','manual-masks')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specify which models and results to include\n",
    "Making changes here may or may not break subsequent as it relies on some of these things being present (e.g., R2, R2_1, and R2_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Models = {\n",
    "    'dhrf' : ['linhrf_cv1_dhrf','linhrf_cv1_dhrf_neggain','csshrf_cv1_dhrf','doghrf_cv1_dhrf'],\n",
    "    'mhrf' : ['linhrf_cv1_mhrf','linhrf_cv1_mhrf_neggain','csshrf_cv1_mhrf','doghrf_cv1_mhrf'],\n",
    "    'names' : ['lin','lin_ng','css','dog'],\n",
    "}\n",
    "Res_type = ['ANG', 'ECC', 'EXPT', 'FWHM', 'GAIN', 'IMAG', 'REAL', 'RFS', 'X', 'Y', 'NAMP', 'SDRATIO']\n",
    "xfm = 'epi2surf'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a volume to vertex mapper for this subject"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vol2Vert={}\n",
    "for s in subs:\n",
    "    Vol2Vert[s] = cortex.get_mapper(s, xfm, 'line_nearest', recache=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a dictionary to collect all FitResults in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this will be the structure of data collection\n",
    "SR={}\n",
    "for s in subs:\n",
    "    SR[s] = {\n",
    "        'subject' : s,\n",
    "        'xfm' : xfm,\n",
    "        'mHRF' : {\n",
    "            'arr' : {},\n",
    "            'vol' : {},\n",
    "        },\n",
    "        'dHRF' : {\n",
    "            'arr' : {},\n",
    "            'vol' : {},\n",
    "        },\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Populate the dictionary with the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the data\n",
    "for s in subs:\n",
    "    for h in ['m','d']:\n",
    "        print('Processing ' + h + 'HRF:', end='')\n",
    "        for M in Models[h + 'hrf']:\n",
    "            # get model idx so we create a shorter variable name\n",
    "            midx = Models[h + 'hrf'].index(M)\n",
    "            print(' ' + Models['names'][midx], end='')\n",
    "\n",
    "            # get the mean R2 map\n",
    "            volpath = os.path.join(FitResPath,'MRI',s.lower(),M,'Sess-' + M + '_meanR2.nii.gz')\n",
    "            volpath1 = os.path.join(FitResPath,'MRI',s.lower(),M,'Sess-' + M + '_R2_1.nii.gz')\n",
    "            volpath2 = os.path.join(FitResPath,'MRI',s.lower(),M,'Sess-' + M + '_R2_2.nii.gz')\n",
    "\n",
    "            # load the results into numpy arrays with nibabel\n",
    "            R2 = np.array(nib.load(volpath).dataobj)\n",
    "            R2_1 = np.array(nib.load(volpath1).dataobj)\n",
    "            R2_2 = np.array(nib.load(volpath2).dataobj)\n",
    "\n",
    "            # convert to pycortex volumes\n",
    "            R2v = cortex.Volume(R2.transpose(2,1,0), s, xfm)\n",
    "            R2_1v = cortex.Volume(R2_1.transpose(2,1,0), s, xfm)\n",
    "            R2_2v = cortex.Volume(R2_2.transpose(2,1,0), s, xfm)\n",
    "\n",
    "            # add info to dictionaries\n",
    "            # numpy arrays\n",
    "            FitRes = {\n",
    "                'R2' : R2,\n",
    "                'R2_1' : R2_1,\n",
    "                'R2_2' : R2_2\n",
    "                }\n",
    "            # pycortex volumes\n",
    "            FitRes_vol = {\n",
    "                'R2' : R2v,\n",
    "                'R2_1' : R2_1v,\n",
    "                'R2_2' : R2_2v\n",
    "                }\n",
    "\n",
    "            # also get othere results\n",
    "            for R in Res_type:\n",
    "                volpath = os.path.join(FitResPath,'MRI',s.lower(),M,'TH_0', R + '_th0.nii.gz')\n",
    "                if os.path.exists(volpath):\n",
    "                    FitRes[R] = np.array(nib.load(volpath).dataobj)\n",
    "                    FitRes_vol[R] = cortex.Volume(FitRes[R].transpose(2,1,0), s, xfm)\n",
    "\n",
    "            # bring it all together\n",
    "            SR[s][h + 'HRF']['arr'][Models['names'][midx]] = FitRes\n",
    "            SR[s][h + 'HRF']['vol'][Models['names'][midx]] = FitRes_vol\n",
    "        print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if False:\n",
    "    subj = 'M02'\n",
    "    # Check in webviewer \n",
    "    cortex.webgl.show(data=SR[subj]['mHRF']['vol']['lin'])\n",
    "    cortex.webgl.show(data=SR[subj]['mHRF']['vol']['lin_ng'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R2_diffarray={}\n",
    "for s in subs:\n",
    "    R2_diffarray[s] = SR[s]['mHRF']['arr']['lin']['R2'] - SR[s]['mHRF']['arr']['lin_ng']['R2']\n",
    "    SR[s]['mHRF']['vol']['lin_ng']['R2_DIFF'] = cortex.Volume(R2_diffarray[s].transpose(2,1,0), s, xfm)\n",
    "    #cortex.webgl.show(data=SR[s]['mHRF']['vol']['lin_ng'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do some checks to see if this worked as expected\n",
    "Once we know it worked we can switch this off again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if False:\n",
    "    subj='M01'\n",
    "    # Inspect the result volumes (should be numpy arrays)\n",
    "    RR = copy.copy(SR[subj]['mHRF']['arr']['lin']['R2']) # copy the R2 values for some model\n",
    "    RR[RR < 5] ='nan' # threshold it to some level\n",
    "    \n",
    "    cortex.webgl.show(data=cortex.Volume(RR.transpose(2,1,0), subj, xfm)) # check to seeit worked"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function to mask the results based on R2 values\n",
    "This function makes a deepcopy of the input dictionary to prevend overwriting it with masked data. Than masks it with the provided R2-threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def r2mask(DataDict,RTH):\n",
    "    DD = copy.deepcopy(DataDict) # copy the original data so it won't get overwritten\n",
    "    # mask the numpy arrays by inserting nan's\n",
    "    for hrf in ['mHRF','dHRF']:\n",
    "        for m in DD[hrf]['arr']:\n",
    "            # mask all available outputs except for R2\n",
    "            for res in DD[hrf]['arr'][m]:\n",
    "                if res is not 'R2':\n",
    "                    with np.errstate(invalid='ignore'):\n",
    "                        DD[hrf]['arr'][m][res][ DD[hrf]['arr'][m]['R2'] < RTH ] = 'nan'\n",
    "                    # also convert to the pycortex volume\n",
    "                    DD[hrf]['vol'][m][res] = cortex.Volume(DD[hrf]['arr'][m][res].transpose(2,1,0), DD['subject'], DD['xfm'])\n",
    "            # mask R2            \n",
    "            with np.errstate(invalid='ignore'):\n",
    "                DD[hrf]['arr'][m]['R2'][ DD[hrf]['arr'][m]['R2'] < RTH ] = 'nan' \n",
    "            DD[hrf]['vol'][m]['R2m'] = cortex.Volume(DD[hrf]['arr'][m]['R2'].transpose(2,1,0), DD['subject'], DD['xfm'])\n",
    "\n",
    "    return DD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check whether the masking function works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if False:\n",
    "    mFR = r2mask(SR[subj],2)\n",
    "    cortex.webgl.show(data=mFR['mHRF']['vol']['lin'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get D99 atlas information "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read atlas labels\n",
    "D99_inFunc = {}\n",
    "for subj in subs:\n",
    "    D99_inFunc[subj]={}\n",
    "    D99_inFunc[subj]['path'] = os.path.join(ManualMaskPath,'sub-' + subj.lower(),'atlas','D99_in_' + subj + '_adj_inFunc.nii')\n",
    "    D99_inFunc[subj]['arr'] = np.array(nib.load(D99_inFunc[subj]['path']).dataobj)\n",
    "    D99_inFunc[subj]['vol'] = cortex.Volume(D99_inFunc[subj]['arr'].transpose(2,1,0), subj, xfm)\n",
    "    D99_inFunc[subj]['labelpath'] = os.path.join(ManualMaskPath,'sub-' + subj.lower(),'atlas','D99_labeltable_reformat.txt')\n",
    "\n",
    "    D99_inFunc[subj]['labels'] = {}\n",
    "    with open(D99_inFunc[subj]['labelpath']) as f:\n",
    "        for line in f:\n",
    "            labelnum, label = line.strip().split(' ',1)\n",
    "            D99_inFunc[subj]['labels'][label.strip()] = int(labelnum)\n",
    "    #print(D99_inFunc['labels'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if False:\n",
    "    subj='M01'\n",
    "    cortex.webgl.show(data=D99_inFunc[subj]['vol'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function that returns ROI names from voxel label number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_roiname(LabelDict,LabelNum):\n",
    "    for name, number in LabelDict.items():\n",
    "        if number == LabelNum:\n",
    "            return name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test if the get_roiname function works\n",
    "if False:\n",
    "    subj='M01'\n",
    "    roi = get_roiname(D99_inFunc[subj]['labels'],147)\n",
    "    print(roi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function that returns data masked by ROI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roimask(DataDict,AtlasDict,rois):\n",
    "    # DataDict is the dictionary of results\n",
    "    # AtlasDict is a dictionary of D99 atlas info\n",
    "    # rois is a list of rois to include\n",
    "    \n",
    "    DD = copy.deepcopy(DataDict) # copy the original data so it won't get overwritten\n",
    "    AA = copy.deepcopy(AtlasDict) # copy the original data so it won't get overwritten\n",
    "    \n",
    "    # Create a mask that includes all rois in the list\n",
    "    mask = np.zeros(AA['arr'].shape)\n",
    "    for r in rois:\n",
    "        tmask = AA['arr']==AA['labels'][r]\n",
    "        mask = mask + tmask\n",
    "    mask[mask > 0]\n",
    "    \n",
    "    # mask the numpy arrays by inserting nan's\n",
    "    for hrf in ['mHRF','dHRF']:\n",
    "        for m in DD[hrf]['arr']:\n",
    "            # mask all available outputs except for R2\n",
    "            for res in DD[hrf]['arr'][m]:\n",
    "                with np.errstate(invalid='ignore'):\n",
    "                    DD[hrf]['arr'][m][res][ mask < 1 ] = 'nan'\n",
    "                # also convert to the pycortex volume\n",
    "                DD[hrf]['vol'][m][res] = cortex.Volume(DD[hrf]['arr'][m][res].transpose(2,1,0), DD['subject'], DD['xfm'])\n",
    "    return DD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test this function\n",
    "if False:\n",
    "    subj='M01'\n",
    "    V1Res = roimask(SR[subj],D99_inFunc[subj],['V1'])\n",
    "    cortex.webgl.show(data=V1Res['mHRF']['vol']['lin'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test masking by ROI and R2\n",
    "if False:\n",
    "    subj='M01'\n",
    "    V1Res = roimask(SR[subj],D99_inFunc[subj],['V1'])\n",
    "    mV1Res = r2mask(V1Res,4)\n",
    "    cortex.webgl.show(data=mV1Res['mHRF']['vol']['lin'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create visualizations for manuscript figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subj = 'M02'\n",
    "hrf = 'mHRF'\n",
    "model = 'css'\n",
    "\n",
    "M1 = r2mask(SR[subj],5)\n",
    "#cortex.webgl.show(data=M1[hrf]['vol'][model])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect the localization of negative model benefits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TH = 5\n",
    "hrf = 'mHRF'\n",
    "\n",
    "DoG={}; PLIN={}; ULIN={}\n",
    "   \n",
    "# NB! Rememeber that masks are 'inverse masks' here, used to set the non-wanted voxels to 'nan'       \n",
    "    \n",
    "for s in subs:  \n",
    "    DoG[s]={}\n",
    "    PLIN[s]={}\n",
    "    ULIN[s]={}\n",
    "    \n",
    "    DoG[s]['R2'] = copy.deepcopy(SR[s][hrf]['arr']['dog']['R2'])\n",
    "    PLIN[s]['R2'] = copy.deepcopy(SR[s][hrf]['arr']['lin']['R2'])\n",
    "    ULIN[s]['R2'] = copy.deepcopy(SR[s][hrf]['arr']['lin_ng']['R2'])\n",
    "\n",
    "    ULIN[s]['R2_above_th'] = copy.copy(ULIN[s]['R2'])\n",
    "    ULIN[s]['R2_sub_th'] = copy.copy(ULIN[s]['R2'])\n",
    "    ULIN[s]['THMASK'] = ULIN[s]['R2'] > TH\n",
    "    ULIN[s]['R2_above_th'][ np.invert(ULIN[s]['THMASK']) ] = 'nan'\n",
    "    ULIN[s]['R2_sub_th'][ ULIN[s]['THMASK'] ] = 'nan'\n",
    "    \n",
    "    PLIN[s]['R2_above_th'] = copy.copy(PLIN[s]['R2'])\n",
    "    PLIN[s]['R2_sub_th'] = copy.copy(PLIN[s]['R2'])\n",
    "    PLIN[s]['THMASK'] = PLIN[s]['R2'] > TH\n",
    "    PLIN[s]['R2_above_th'][ np.invert(PLIN[s]['THMASK']) ] = 'nan'\n",
    "    PLIN[s]['R2_sub_th'][ PLIN[s]['THMASK'] ] = 'nan'\n",
    "\n",
    "    DoG[s]['R2_above_th'] = copy.copy(DoG[s]['R2'])\n",
    "    DoG[s]['R2_sub_th'] = copy.copy(DoG[s]['R2'])\n",
    "    DoG[s]['THMASK'] = DoG[s]['R2'] > TH\n",
    "    DoG[s]['R2_above_th'][ np.invert(DoG[s]['THMASK']) ] = 'nan'\n",
    "    DoG[s]['R2_sub_th'][ DoG[s]['THMASK'] ] = 'nan'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect U-LIN gain and DoG nAMP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "UD={}; DN={};\n",
    "for s in subs:\n",
    "    UD[s]={}\n",
    "    UD[s]['Gain'] = SR[s][hrf]['arr']['lin_ng']['GAIN']\n",
    "    UD[s]['ULIN_R2'] = ULIN[s]['R2_above_th']\n",
    "    UD[s]['Gain_masked'] = copy.copy(UD[s]['Gain'])\n",
    "    UD[s]['Gain_masked'][np.invert(ULIN[s]['THMASK'])] = 'nan'\n",
    "    \n",
    "    UD[s]['nAmp'] = SR[s][hrf]['arr']['dog']['NAMP']\n",
    "    UD[s]['DoG_R2'] = DoG[s]['R2_above_th']\n",
    "    UD[s]['nAmp_masked'] = copy.copy(UD[s]['nAmp'])\n",
    "    UD[s]['nAmp_masked'][np.invert(DoG[s]['THMASK'])] = 'nan'  \n",
    "    \n",
    "    \n",
    "    for a in UD[s]:\n",
    "        print(a)\n",
    "        DN[s + '_' + a] = cortex.Volume(UD[s][a].transpose(2,1,0), s, xfm)\n",
    "        \n",
    "cortex.webgl.show(data=DN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## U-LIN vs P-LIN    \n",
    "Here we check the voxels for which both U-LIN and P-LIN fit decently"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# P-LIN & U-LIN R2 > TH\n",
    "# Map dR2 U-LIN minus P-LIN\n",
    "hrf = 'mHRF'\n",
    "UP={}; DN={};\n",
    "for s in subs:\n",
    "    UP[s]={}\n",
    "    UP[s]['dR2_th'] = ULIN[s]['R2_above_th'] - PLIN[s]['R2_above_th']    \n",
    "    UP[s]['nanMASK'] = np.isnan(UP[s]['dR2_th'])\n",
    "    \n",
    "    UP[s]['Ecc'] = SR[s][hrf]['arr']['lin']['ECC'] \n",
    "    UP[s]['Ecc_masked'] = copy.copy(UP[s]['Ecc'])\n",
    "    UP[s]['Ecc_masked'][UP[s]['nanMASK']] = 'nan'\n",
    "    \n",
    "    UP[s]['Ecc_ng'] = SR[s][hrf]['arr']['lin_ng']['ECC'] \n",
    "    UP[s]['Ecc_ng_masked'] = copy.copy(UP[s]['Ecc_ng'])\n",
    "    UP[s]['Ecc_ng_masked'][UP[s]['nanMASK']] = 'nan'\n",
    "    \n",
    "    for a in UP[s]:\n",
    "        print(a)\n",
    "        DN[s + '_' + a] = cortex.Volume(UP[s][a].transpose(2,1,0), s, xfm)\n",
    "        \n",
    "cortex.webgl.show(data=DN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## U-LIN vs P-LIN     \n",
    "Here we check the voxels for which U-LIN performs well, but P-LIN does not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# U-LIN R2 > TH & P-LIN < TH\n",
    "# Map dR2 U-LIN minus P-LIN\n",
    "UP={}; DN={};\n",
    "for s in subs:\n",
    "    UP[s]={}\n",
    "    UP[s]['dR2_th'] = ULIN[s]['R2_above_th'] - PLIN[s]['R2_sub_th']\n",
    "    #UP[s]['ULIN-MASK'] = np.invert( ULIN[s]['THMASK'] )\n",
    "    #UP[s]['PLIN-MASK'] = PLIN[s]['THMASK']    \n",
    "    UP[s]['nanMASK'] = np.isnan(UP[s]['dR2_th'])\n",
    "    \n",
    "    UP[s]['Ecc'] = SR[s][hrf]['arr']['lin']['ECC'] \n",
    "    UP[s]['Ecc_masked'] = copy.copy(UP[s]['Ecc'])\n",
    "    UP[s]['Ecc_masked'][UP[s]['nanMASK']] = 'nan'\n",
    "\n",
    "    UP[s]['R2_ULIN'] = ULIN[s]['R2']\n",
    "    UP[s]['R2_PLIN'] = PLIN[s]['R2']\n",
    "    \n",
    "    for a in UP[s]:\n",
    "        print(a)\n",
    "        DN[s + '_' + a] = cortex.Volume(UP[s][a].transpose(2,1,0), s, xfm)\n",
    "        \n",
    "cortex.webgl.show(data=DN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DoG vs P-LIN    \n",
    "Here we check the voxels for which DoG outperforms P-LIN, while both perform ok."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DoG R2 > TH & P-LIN > TH\n",
    "# Map dR2 DoG minus P-LIN\n",
    "UP={}; DN={};\n",
    "for s in subs:\n",
    "    UP[s]={}\n",
    "    UP[s]['dR2_th'] = DoG[s]['R2_above_th'] - PLIN[s]['R2_above_th'] \n",
    "    UP[s]['nanMASK'] = np.isnan(UP[s]['dR2_th'])\n",
    "    \n",
    "    UP[s]['Ecc'] = SR[s][hrf]['arr']['lin']['ECC'] \n",
    "    UP[s]['Ecc_masked'] = copy.copy(UP[s]['Ecc'])\n",
    "    UP[s]['Ecc_masked'][UP[s]['nanMASK']] = 'nan'\n",
    "        \n",
    "    UP[s]['R2_DoG'] = DoG[s]['R2']\n",
    "    UP[s]['R2_PLIN'] = PLIN[s]['R2']\n",
    "    \n",
    "    for a in UP[s]:\n",
    "        print(a)\n",
    "        DN[s + '_' + a] = cortex.Volume(UP[s][a].transpose(2,1,0), s, xfm)\n",
    "        \n",
    "cortex.webgl.show(data=DN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}