3_create_figures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figures\n",
    "\n",
    "This notebook creates the following figures:\n",
    "\n",
    "* visual field coverage plots\n",
    "* polar angle histograms\n",
    "* eccentricity distributions\n",
    "* eccentricity-size relations\n",
    "* timeseries plot\n",
    "* across-split correlation of parameters\n",
    "* tSNR - R2 correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# imports\n",
    "from __future__ import division\n",
    "import os\n",
    "import nibabel as nb\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import h5py\n",
    "import matplotlib.pyplot as pl\n",
    "import seaborn as sn\n",
    "sn.set_style('ticks')\n",
    "import glob\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.stats.api as sms\n",
    "from tqdm import tqdm\n",
    "\n",
    "\n",
    "import pandas as pd\n",
    "import shutil\n",
    "import pycircstat as pc\n",
    "\n",
    "import matplotlib as mpl\n",
    "%matplotlib inline\n",
    "mpl.rc_file_defaults()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### define directories:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# set this to your data directory:\n",
    "data_dir = '/home/shared/2018/visual/hcp_cerebellum/'\n",
    "repo_dir = '/home/vanes/git/hcp_cerebellum/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### subdirectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "resource_dir = os.path.join(repo_dir,'resources')\n",
    "retmap_dir = os.path.join(resource_dir,'volume_masks')\n",
    "\n",
    "# setup figure directory\n",
    "fig_dir = os.path.join(data_dir,'figs')\n",
    "if not os.path.isdir(fig_dir): os.mkdir(fig_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# these are the dimensions of the niftis\n",
    "dims = {\n",
    "    'ang':0,\n",
    "    'ecc':1,\n",
    "    'gain':2,\n",
    "    'meanvol':3,\n",
    "    'r2':4,\n",
    "    'rfsize':5\n",
    "}   \n",
    "\n",
    "stim_radius = 8\n",
    "avg_masktype = 'r2_spill_fix'\n",
    "ind_masktype = 'r2_fix_roi'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### define rois"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# these refer to the indices in cerebellum_retmaps.nii\n",
    "mask_names ={\n",
    "'left_mOMV':7,\n",
    "'right_mOMV':9,\n",
    "'left_lOMV':8,\n",
    "'right_lOMV':10,\n",
    "    \n",
    "'left_VIIb':5,\n",
    "'right_VIIb':6,\n",
    "    \n",
    "'left_mVIIIb':3,\n",
    "'right_mVIIIb':2,\n",
    "'left_lVIIIb':4,\n",
    "'right_lVIIIb':1,\n",
    "}\n",
    "\n",
    "# this is how the lateral/medial and left right hemisphere indices combine:\n",
    "roi_combs = {\n",
    "    'left_OMV':[7,8],\n",
    "    'right_OMV':[9,10],\n",
    "    'left_VIIb':[5],\n",
    "    'right_VIIb':[6],\n",
    "    'left_VIIIb':[3,4],\n",
    "    'right_VIIIb':[1,2],\n",
    "\n",
    "    'OMV':[7,8,9,10],\n",
    "    'lOMV':[8,10],\n",
    "    'mOMV':[7,9],\n",
    "    'VIIb':[5,6],\n",
    "    'VIIIb':[1,2,3,4],\n",
    "    'lVIIIb':[1,4],\n",
    "    'mVIIIb':[2,3],\n",
    "    \n",
    "    'all_rois':[1,2,3,4,5,6,7,8,9,10]\n",
    "}\n",
    "\n",
    "roi_order = ['left_mOMV','right_mOMV','left_lOMV','right_lOMV','left_VIIb','right_VIIb','left_mVIIIb','right_mVIIIb','left_lVIIIb','right_lVIIIb']\n",
    "roi_comb_order = ['mOMV','lOMV','VIIb','mVIIIb','lVIIIb']           "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get the data\n",
    "sjs = range(181) + ['181','181_1','181_2','182','182_1','182_2','183','183_1','183_2']\n",
    "data = {}\n",
    "for sj in sjs:\n",
    "    if sj == '183':\n",
    "        data[sj] =nb.load(os.path.join(data_dir,'masked_niftis',avg_masktype,'prfresults_subject_rank_%s.nii.gz'%sj)).get_data() \n",
    "    else:\n",
    "        data[sj] =nb.load(os.path.join(data_dir,'masked_niftis',ind_masktype,'prfresults_subject_rank_%s.nii.gz'%sj)).get_data() \n",
    "    \n",
    "cer_retmaps = nb.load(os.path.join(retmap_dir,'cerebellum_retmaps.nii')).get_data()\n",
    "cer_mask = nb.load(os.path.join(retmap_dir,'cmask.nii')).get_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### function for determining roi color"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_roi_color(mask):\n",
    "    \n",
    "    # determine base roi colors\n",
    "    if mask in ['left_OMV','right_OMV','OMV']:\n",
    "        c = '#6590CB'\n",
    "    elif mask in ['left_VIIIb','right_VIIIb','VIIIb']:\n",
    "        c = '#E55D5C'\n",
    "    elif ('VIIb' in mask) + (mask=='VIIb'):\n",
    "        c = '#E79F2A'\n",
    "        \n",
    "    # roi colors for lateral / medial separation:\n",
    "    elif mask in ['left_lOMV','right_lOMV','lOMV']:\n",
    "        c = '#6590CB'\n",
    "    elif mask in ['left_mOMV','right_mOMV','mOMV']:\n",
    "        c = '#0D5B99'\n",
    "    elif mask in ['left_mVIIIb','right_mVIIIb','mVIIIb']:\n",
    "        c = '#C03142'\n",
    "    elif mask in ['left_lVIIIb','right_lVIIIb','lVIIIb']:\n",
    "        c = '#E57F80'\n",
    "        \n",
    "    # for all rois and masks\n",
    "    elif mask in ['all_rois']:\n",
    "        c = 'k'\n",
    "    elif mask in ['all_voxels']:\n",
    "        c = 'gray'\n",
    "        \n",
    "    return c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualize pRF center + size distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for sj in ['181','182','183','183_1','183_2']:\n",
    "    avgdata = data[sj]\n",
    "    f = pl.figure(figsize=(2,4.5))\n",
    "    for mi,mask in enumerate(roi_order):\n",
    "\n",
    "        s = f.add_subplot(5,2,mi+1,aspect='equal')\n",
    "\n",
    "        # determine this roi mask\n",
    "        roimask = (cer_retmaps==(mask_names[mask]))\n",
    "\n",
    "        # convert polar to cartesian:\n",
    "        xs = np.cos(np.radians(np.ravel(avgdata[roimask,dims['ang']]))) * np.ravel(avgdata[roimask,dims['ecc']])\n",
    "        ys = np.sin(np.radians(np.ravel(avgdata[roimask,dims['ang']]))) * np.ravel(avgdata[roimask,dims['ecc']])\n",
    "        eccs = np.ravel(avgdata[roimask,dims['ecc']])\n",
    "\n",
    "        # get sizes\n",
    "        sizes = np.ravel(avgdata[roimask,dims['rfsize']])\n",
    "        \n",
    "        # generate mask\n",
    "        v = ~np.isnan(sizes)\n",
    "        v *= (eccs < stim_radius)\n",
    "        \n",
    "        # apply mask\n",
    "        xs = xs[v]\n",
    "        ys = ys[v]\n",
    "        sizes = sizes[v]\n",
    "\n",
    "        # determine roi color\n",
    "        c = get_roi_color(mask)\n",
    "\n",
    "        # now draw a circle at rf size\n",
    "        for x,y,sigma in zip(xs,ys,sizes):\n",
    "            s.add_artist(pl.Circle((x,y),sigma, color=c,fill=False,alpha=0.25))\n",
    "\n",
    "        # draw crosshair\n",
    "        pl.axhline(0,lw=0.5,color='k')\n",
    "        pl.axvline(0,lw=0.5,color='k')\n",
    "\n",
    "        # draw dot at center \n",
    "        pl.plot(xs,ys,'o',color=c,ms=3,mec='w',mew=1,alpha=1)\n",
    "        pl.xlim(-stim_radius,stim_radius)\n",
    "        pl.ylim(-stim_radius,stim_radius)   \n",
    "\n",
    "        # only draw the axis labels for VIIIb:\n",
    "        if mask =='left_lVIIIb':\n",
    "            pl.xticks([-stim_radius,0,stim_radius],['-%d'%stim_radius,0,'%.d'%stim_radius])\n",
    "            pl.yticks([-stim_radius,0,stim_radius],['-%d'%stim_radius,0,'%.d'%stim_radius])\n",
    "            sn.despine(offset=2)\n",
    "        elif mask =='right_lVIIIb':\n",
    "            pl.xticks([-stim_radius,0,stim_radius],['-%.d'%stim_radius,0,'%.d'%stim_radius])\n",
    "            pl.yticks([])\n",
    "            sn.despine(offset=2,left=True)\n",
    "        else:\n",
    "            pl.axis('off')\n",
    "            pl.xticks([])\n",
    "            pl.yticks([])                \n",
    "\n",
    "    pl.tight_layout()\n",
    "    f.savefig(os.path.join(fig_dir,'prf_scatter_avg_subject_%s.pdf'%sj))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### create eccentricity-size plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for sj in ['181','182','183','183_1','183_2']:\n",
    "    avgdata = data[sj]\n",
    "    # eccen-size plot\n",
    "    this_roi_comb_order = ['mVIIIb','lVIIIb','VIIb','mOMV','lOMV']\n",
    "\n",
    "    f = pl.figure(figsize=(1.25,1.25))\n",
    "    s = f.add_subplot(111)\n",
    "\n",
    "    all_max_eccs = []\n",
    "    all_max_sizes=[]\n",
    "    for mi,mask in enumerate(this_roi_comb_order):\n",
    "\n",
    "        # determine this roi mask\n",
    "        roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "        for subroi in roi_combs[mask]:\n",
    "            roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "        # determine roi color\n",
    "        c = get_roi_color(mask)  \n",
    "\n",
    "        sizes = np.ravel(avgdata[roimask,dims['rfsize']])\n",
    "        eccs = np.ravel(avgdata[roimask,dims['ecc']])\n",
    "\n",
    "        maxecc = np.min([np.nanmax(eccs),stim_radius])\n",
    "    \n",
    "        # generate mask\n",
    "        v = (eccs < maxecc)\n",
    "        v*= ~np.isnan(sizes)\n",
    "\n",
    "        # do linear regression with statsmodels\n",
    "        y = sizes[v]\n",
    "        X = eccs[v]\n",
    "        X = sm.add_constant(X)     \n",
    "\n",
    "        mod = sm.OLS(y, X)\n",
    "        res = mod.fit()\n",
    "        ci_intercept,ci_slope = res.conf_int(0.05)\n",
    "        mean_intercept,mean_slope = res.params\n",
    "\n",
    "        fit = np.polyval([mean_slope,mean_intercept],np.linspace(0,maxecc,15))\n",
    "        fit_high = np.polyval([ci_slope[0],ci_intercept[1]],np.linspace(0,maxecc,15))\n",
    "        fit_low = np.polyval([ci_slope[1],ci_intercept[0]],np.linspace(0,maxecc,15))\n",
    "\n",
    "        pl.plot(np.linspace(0,maxecc,15),fit,color=c,lw=2)\n",
    "        pl.fill_between(np.linspace(0,maxecc,15),fit_low,fit_high,alpha=0.25,color=c)\n",
    "\n",
    "    pl.xlim(0,stim_radius)  \n",
    "    pl.ylim(0,10)\n",
    "    pl.xticks([0,stim_radius])\n",
    "    pl.yticks([0,10])\n",
    "\n",
    "    sn.despine(offset=2)\n",
    "\n",
    "    pl.tight_layout()\n",
    "    f.savefig(os.path.join(fig_dir,'prf_eccsize_avg_%s.pdf'%sj))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### create eccentricity histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for sj in ['181','182','183','183_1','183_2']:\n",
    "    avgdata = data[sj]\n",
    "    f = pl.figure(figsize=(1.25,2))\n",
    "    s = f.add_subplot(111)\n",
    "\n",
    "    dc = pd.DataFrame()\n",
    "    colors = []\n",
    "\n",
    "    for mi,mask in enumerate(roi_comb_order):\n",
    "\n",
    "        # determine this roi mask\n",
    "        roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "        for subroi in roi_combs[mask]:\n",
    "            roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "        # get data\n",
    "        eccs = np.ravel(avgdata[roimask,dims['ecc']])\n",
    "        \n",
    "        v = ~np.isnan(eccs)\n",
    "        \n",
    "        eccs = eccs[v]\n",
    "\n",
    "        # get roi color and append to all\n",
    "        c = get_roi_color(mask)\n",
    "        colors.append(c)\n",
    "\n",
    "        # get eccentricities, append with nans (otherwise df is shortened to roi with least voxels)\n",
    "        dc[mask] = np.hstack([eccs,np.ones(100000-len(eccs))*np.nan])\n",
    "\n",
    "    sn.violinplot(data=dc,palette=colors,scale='width',orient='h',linewidth=0)\n",
    "\n",
    "    pl.xlim(0,stim_radius)\n",
    "    pl.xticks([0,stim_radius])\n",
    "    sn.despine(offset=2)\n",
    "    pl.xlabel('eccen (dva)')\n",
    "\n",
    "    pl.tight_layout(pad=0)\n",
    "    f.savefig(os.path.join(fig_dir,'ecc_histograms_avg_%s.pdf'%sj))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### polar histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for sj in ['181','182','183','183_1','183_2']:\n",
    "    avgdata = data[sj]\n",
    "    for hemi in ['left','right']:\n",
    "        f = pl.figure(figsize=(0.75,4.5))\n",
    "\n",
    "        for mi,roi in enumerate(['mOMV','lOMV','VIIb','mVIIIb','lVIIIb']):\n",
    "            mask = hemi+'_'+roi\n",
    "            s = f.add_subplot(5,1,mi+1,projection='polar')\n",
    "\n",
    "            # determine this roi mask\n",
    "            roimask = (cer_retmaps==(mask_names[mask]))\n",
    "\n",
    "            # determine color for this ROI\n",
    "            c = get_roi_color(mask)\n",
    "\n",
    "            # get data (convert to radians for polar plot)\n",
    "            angles = np.radians(np.ravel(avgdata[roimask,dims['ang']]))\n",
    "            eccs = np.ravel(avgdata[roimask,dims['ecc']])\n",
    "\n",
    "            v = ~np.isnan(angles)\n",
    "            v *= (eccs < stim_radius)\n",
    "            angles = angles[v]\n",
    "            \n",
    "            # compute bins (let's do for 12 bins)\n",
    "            bins_number = 12\n",
    "            bins = np.linspace(0.0, 2 * np.pi, bins_number + 1)\n",
    "            n, _, _ = pl.hist(angles, bins)\n",
    "            width = (2 * np.pi) / bins_number\n",
    "\n",
    "            # plot polar bars\n",
    "            s.cla()        \n",
    "            pl.bar(bins[:bins_number], n, width=width, bottom=0.0,ec=c,color=c)\n",
    "\n",
    "            # plot settings\n",
    "            pl.xticks([])\n",
    "            pl.yticks([])\n",
    "            pl.ylim(0,np.max(n))\n",
    "\n",
    "        pl.tight_layout()\n",
    "        f.savefig(os.path.join(fig_dir,'polar_histograms_avg_%s_%s.pdf'%(hemi,sj)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### individual subject eccen-size relations - selected subjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# eccen-size plot\n",
    "for sgi,sjg in zip([range(5),range(88,93),range(176,181)],['top','mid','bottom']):\n",
    "    \n",
    "    roi_comb_order = ['mOMV','lOMV','VIIb','mVIIIb','lVIIIb']               \n",
    "\n",
    "\n",
    "    f = pl.figure(figsize=(0.75*5,0.75))\n",
    "    for mi,mask in enumerate(roi_comb_order):\n",
    "        s = f.add_subplot(1,len(roi_comb_order),mi+1)\n",
    "\n",
    "        # determine this roi mask\n",
    "        roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "        for subroi in roi_combs[mask]:\n",
    "            roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "        # determine ecc limit and roi color\n",
    "        c = get_roi_color(mask)  \n",
    "\n",
    "        for sj in sgi:\n",
    "            \n",
    "            this_sj_data = nb.load(os.path.join(data_dir,'masked_niftis',ind_masktype,'prfresults_subject_rank_%d.nii.gz'%sj)).get_data()\n",
    "\n",
    "            try:                \n",
    "                sizes = np.ravel(this_sj_data[roimask,dims['rfsize']])\n",
    "                eccs = np.ravel(this_sj_data[roimask,dims['ecc']])                \n",
    "                \n",
    "                maxecc = np.min([np.nanmax(eccs),stim_radius])\n",
    "\n",
    "                # determine valid voxels\n",
    "                v = (eccs < maxecc)\n",
    "                v *= ~np.isnan(eccs)\n",
    "\n",
    "                # do linear regression with statsmodels\n",
    "                y = sizes[v]\n",
    "                X = eccs[v]\n",
    "                X = sm.add_constant(X)     \n",
    "\n",
    "                mod = sm.OLS(y, X)\n",
    "                res = mod.fit()\n",
    "                ci_intercept,ci_slope = res.conf_int(0.05)\n",
    "                mean_intercept,mean_slope = res.params\n",
    "\n",
    "\n",
    "                fit = np.polyval([mean_slope,mean_intercept],np.linspace(0,maxecc,15))\n",
    "                fit_high = np.polyval([ci_slope[0],ci_intercept[1]],np.linspace(0,maxecc,15))\n",
    "                fit_low = np.polyval([ci_slope[1],ci_intercept[0]],np.linspace(0,maxecc,15))\n",
    "\n",
    "                pl.plot(np.linspace(0,maxecc,15),fit,color=c,lw=1)\n",
    "                pl.fill_between(np.linspace(0,maxecc,15),fit_low,fit_high,alpha=0.25,color=c)\n",
    "                \n",
    "            except:\n",
    "#                 print 'could not fit %s ecc-size relation for %s subject %d (too little data)'%(mask,sjg,sj)\n",
    "                pass\n",
    "\n",
    "        pl.xlim(0,stim_radius)  \n",
    "        pl.ylim(0,10)\n",
    "        pl.xticks([0,stim_radius])\n",
    "        pl.yticks([0,10])\n",
    "\n",
    "        sn.despine(offset=2)\n",
    "\n",
    "    pl.tight_layout()\n",
    "    f.savefig(os.path.join(fig_dir,'prf_eccsize_all_%s.pdf'%(sjg)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### individual subject eccen-size relations - all subjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# eccen-size plot\n",
    "this_roi_comb_order = ['mVIIIb','lVIIIb','VIIb','mOMV','lOMV']\n",
    "min_vox_ns = [5,10,15]\n",
    "for min_vox_n in min_vox_ns:\n",
    "    f = pl.figure(figsize=(1.25,1.25))\n",
    "\n",
    "    for mi,mask in enumerate(this_roi_comb_order):\n",
    "    #     s = f.add_subplot(1,len(roi_comb_order),mi+1)\n",
    "\n",
    "        s = f.add_subplot(1,1,1)\n",
    "\n",
    "        # determine this roi mask\n",
    "        roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "        for subroi in roi_combs[mask]:\n",
    "            roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "        # determine ecc limit and roi color\n",
    "        c = get_roi_color(mask)  \n",
    "\n",
    "        all_slopes = []\n",
    "        all_intercepts = []\n",
    "        maxeccs = []\n",
    "        for sj in range(181):\n",
    "\n",
    "            this_sj_data = nb.load(os.path.join(data_dir,'masked_niftis',ind_masktype,'prfresults_subject_rank_%d.nii.gz'%sj)).get_data()\n",
    "\n",
    "            sizes = np.ravel(this_sj_data[roimask,dims['rfsize']])\n",
    "            eccs = np.ravel(this_sj_data[roimask,dims['ecc']])                \n",
    "\n",
    "            # mask voxels\n",
    "            maxecc = np.min([np.nanmax(eccs),stim_radius])\n",
    "            v = (eccs < maxecc)        \n",
    "            v = ~np.isnan(sizes)\n",
    "\n",
    "            # only go on if there are sufficient voxels in this subject\n",
    "            if v.sum() > min_vox_n:\n",
    "\n",
    "                maxeccs.append(maxecc)\n",
    "\n",
    "                # do linear regression with statsmodels\n",
    "                y = sizes[v]\n",
    "                X = eccs[v]\n",
    "                X = sm.add_constant(X)     \n",
    "\n",
    "                mod = sm.OLS(y, X)\n",
    "                res = mod.fit()\n",
    "                intercept,slope = res.params\n",
    "\n",
    "                all_slopes.append(slope)\n",
    "                all_intercepts.append(intercept)\n",
    "\n",
    "        maxecc = np.mean(maxeccs)\n",
    "\n",
    "        print '%s, n subs: %d, maxecc: %.2f'%(mask,len(all_slopes),maxecc)\n",
    "\n",
    "        # compute mean and ci across valid subjects\n",
    "        mean_slope = np.mean(all_slopes)\n",
    "        mean_intercept = np.mean(all_intercepts)\n",
    "\n",
    "        ci_slope = sms.DescrStatsW(data=np.array(all_slopes)).tconfint_mean()\n",
    "        ci_intercept = sms.DescrStatsW(data=np.array(all_intercepts)).tconfint_mean()\n",
    "    \n",
    "        print 'slope: %.2f, %.2f-%.2f'%(mean_slope,ci_slope[0],ci_slope[1])\n",
    "        print 'intercept: %.2f, %.2f-%.2f'%(mean_intercept,ci_intercept[0],ci_intercept[1])\n",
    "\n",
    "        fit = np.polyval([mean_slope,mean_intercept],np.linspace(0,maxecc,15))\n",
    "        fit_high = np.polyval([ci_slope[0],ci_intercept[1]],np.linspace(0,maxecc,15))\n",
    "        fit_low = np.polyval([ci_slope[1],ci_intercept[0]],np.linspace(0,maxecc,15))\n",
    "\n",
    "        pl.plot(np.linspace(0,maxecc,15),fit,color=c,lw=1)\n",
    "        pl.fill_between(np.linspace(0,maxecc,15),fit_low,fit_high,alpha=0.25,color=c)\n",
    "\n",
    "    pl.xlim(0,stim_radius)  \n",
    "    pl.ylim(0,10)\n",
    "    pl.xticks([0,stim_radius])\n",
    "    pl.yticks([0,10])\n",
    "\n",
    "    sn.despine(offset=2)\n",
    "\n",
    "    pl.tight_layout()\n",
    "    f.savefig(os.path.join(fig_dir,'prf_eccsize_all_subjects_min_voxn_%d.pdf'%min_vox_n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# compare parameters across split half data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def pearsonr_ci(x,y,alpha=0.05):\n",
    "    ''' calculate Pearson correlation along with the confidence interval using scipy and numpy\n",
    "    Parameters\n",
    "    ----------\n",
    "    x, y : iterable object such as a list or np.array\n",
    "      Input for correlation calculation\n",
    "    alpha : float\n",
    "      Significance level. 0.05 by default\n",
    "    Returns\n",
    "    -------\n",
    "    r : float\n",
    "      Pearson's correlation coefficient\n",
    "    pval : float\n",
    "      The corresponding p value\n",
    "    lo, hi : float\n",
    "      The lower and upper bound of confidence intervals\n",
    "      \n",
    "    Taken from: https://zhiyzuo.github.io/Pearson-Correlation-CI-in-Python/\n",
    "    '''\n",
    "    \n",
    "    r, p = sp.stats.pearsonr(x,y)\n",
    "    r_z = np.arctanh(r)\n",
    "    se = 1/np.sqrt(x.size-3)\n",
    "    z = sp.stats.norm.ppf(1-alpha/2)\n",
    "    lo_z, hi_z = r_z-z*se, r_z+z*se\n",
    "    lo, hi = np.tanh((lo_z, hi_z))\n",
    "    return r, p, lo, hi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from statsmodels.stats.weightstats import DescrStatsW\n",
    "\n",
    "def p_val_from_bootstrap_dist(distribution,test_value=0,two_tailed=True):\n",
    "\n",
    "    # see which part of the distribution falls below / above test value:\n",
    "    proportion_smaller_than_test_value = np.sum(np.array(distribution) < test_value) / len(distribution)\n",
    "    proportion_larger_than_test_value = np.sum(np.array(distribution) > test_value) / len(distribution)\n",
    "\n",
    "    # take minimum value as p-val:\n",
    "    p = np.min([proportion_smaller_than_test_value,proportion_larger_than_test_value])\n",
    "\n",
    "    # this yields a one-tailed test, so multiply by 2 if we want a two-tailed p-val:\n",
    "    if two_tailed:\n",
    "        p*=2\n",
    "\n",
    "    return p\n",
    "\n",
    "def get_ci(distribution,ci_factor):\n",
    "\n",
    "    # convert ci factor to percentile\n",
    "    perc_low = (stats.norm.sf(ci_factor))*100.\n",
    "    perc_high = (1 - (stats.norm.sf(ci_factor)))*100.\n",
    "    ci_low = np.percentile(distribution,perc_low)\n",
    "    ci_high = np.percentile(distribution,perc_high)\n",
    "\n",
    "    return [ci_low,ci_high]\n",
    "    \n",
    "def bootstrap_correlation(x_data,y_data,test_value=0,ci_factor = 1.96,two_tailed=True,reps=int(1e3),circ =False):\n",
    "\n",
    "    N = len(x_data)\n",
    "\n",
    "    # get random ints for random indices\n",
    "    permute_indices = np.random.randint(0, len(x_data), size = (len(x_data), int(reps))).T\n",
    "\n",
    "    bootstrap_distr=[]\n",
    "    bootstrap_distr_z=[]\n",
    "    # loop over permutes\n",
    "    for fold in permute_indices:\n",
    "        if circ:\n",
    "            r = pc.corrcc(x_data[fold],y_data[fold])\n",
    "        else:\n",
    "            r = DescrStatsW(data=np.vstack([x_data[fold],y_data[fold]]).T).corrcoef[0,1]\n",
    "        z = np.arctanh(r)\n",
    "        bootstrap_distr.append(r)\n",
    "        bootstrap_distr_z.append(z)\n",
    "\n",
    "    # compute central corr on all data\n",
    "    if circ:\n",
    "        corr = pc.corrcc(x_data,y_data)\n",
    "    else:\n",
    "        corr = DescrStatsW(data=np.vstack([x_data,y_data]).T).corrcoef[0,1]\n",
    "    # calculate p-val\n",
    "\n",
    "    p = p_val_from_bootstrap_dist(bootstrap_distr_z,test_value,two_tailed)\n",
    "\n",
    "\n",
    "    # return standard deviation of bootstrap distro as CI\n",
    "    # corr_ci = np.std(bootstrap_distr)*ci_factor\n",
    "    corr_ci = get_ci(bootstrap_distr,ci_factor)\n",
    "\n",
    "    return corr, p, corr_ci[0],corr_ci[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for comp,name in zip([['183_1','183_2'],['181','182']],['across_runs','across_subjects']):\n",
    "    pl.close('all')\n",
    "\n",
    "    for param in ['ecc','ang','rfsize']:#,'x','y']:\n",
    "        \n",
    "        f = pl.figure(figsize=(1,1))\n",
    "        s = f.add_subplot(111)\n",
    "        pl.axhline(0,color='k')\n",
    "\n",
    "\n",
    "        # eccen-size plot\n",
    "#         this_roi_comb_order = ['mVIIIb','lVIIIb','VIIb','mOMV','lOMV']\n",
    "\n",
    "#         this_roi_comb_order = ['OMV','VIIb','VIIIb']#VIIIb','lVIIIb','VIIb','mOMV','lOMV']\n",
    "        this_roi_comb_order = ['OMV','VIIb','VIIIb']#VIIIb','lVIIIb','VIIb','mOMV','lOMV']\n",
    "\n",
    "\n",
    "        for mi,mask_overall in enumerate(this_roi_comb_order):\n",
    "\n",
    "            for hi,hemi in enumerate(['','_left','_right']):\n",
    "                \n",
    "                mask = hemi + mask_overall\n",
    "\n",
    "                # determine this roi mask\n",
    "                roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "                for subroi in roi_combs[mask]:\n",
    "                    roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "                # determine roi color\n",
    "                c = get_roi_color(mask)  \n",
    "\n",
    "                rfsize1 = data[comp[0]][roimask,dims['rfsize']]\n",
    "                ecc1 = data[comp[0]][roimask,dims['ecc']]\n",
    "                ang1 = data[comp[0]][roimask,dims['ang']]\n",
    "                xs1 = np.cos(np.radians(np.ravel(ang1))) * np.ravel(ecc1)\n",
    "                ys1 = np.sin(np.radians(np.ravel(ang1))) * np.ravel(ecc1)\n",
    "\n",
    "                rfsize2 = data[comp[1]][roimask,dims['rfsize']]\n",
    "                ecc2= data[comp[1]][roimask,dims['ecc']]\n",
    "                ang2 = data[comp[1]][roimask,dims['ang']]\n",
    "\n",
    "                xs2 = np.cos(np.radians(np.ravel(ang2))) * np.ravel(ecc2)\n",
    "                ys2 = np.sin(np.radians(np.ravel(ang2))) * np.ravel(ecc2)\n",
    "\n",
    "\n",
    "                if param == 'x':\n",
    "                    data1 = xs1\n",
    "                    data2 = xs2\n",
    "                    plotrange = [-stim_radius,stim_radius]\n",
    "                elif param == 'y':\n",
    "                    data1 = ys1\n",
    "                    data2 = ys2\n",
    "                    plotrange = [-stim_radius,stim_radius]\n",
    "\n",
    "                elif param == 'ecc':\n",
    "                    data1 = ecc1\n",
    "                    data2 = ecc2\n",
    "                    plotrange = [0,stim_radius]\n",
    "\n",
    "                elif param == 'ang':\n",
    "                    data1 = ang1\n",
    "                    data2 = ang2\n",
    "                    plotrange = [0,360]\n",
    "\n",
    "                elif param == 'rfsize':\n",
    "                    data1 = rfsize1\n",
    "                    data2 = rfsize2\n",
    "                    plotrange = [0,20]\n",
    "\n",
    "\n",
    "                v = (~np.isnan(data1))*(~np.isnan(data2))\n",
    "\n",
    "                data1 = data1[v]\n",
    "                data2 = data2[v]\n",
    "\n",
    "                if param == 'ang':\n",
    "    #                 r,p,l,u = pearsonr_ci(data1,data2)\n",
    "                    r,p,l,u = bootstrap_correlation(np.radians(data1),np.radians(data2),circ=True)\n",
    "                    print(mask,r,p)\n",
    "                else:\n",
    "                    r,p,l,u = bootstrap_correlation(data1,data2)\n",
    "    #                 r,(l,u)= pc.corrcc(np.radians(data1),np.radians(data2),ci=0.95,bootstrap_iter=int(1e4))\n",
    "    #             pl.plot([mi-0.1,mi+0.1],[r,r],c=c,lw=3)\n",
    "\n",
    "                if np.arctanh(l) > 0:\n",
    "                    pl.plot(mi+(0.25*hi),r,'o',ms=2,mec=c,c=c,mew=2)\n",
    "                else:\n",
    "                    pl.plot(mi+(0.25*hi),r,'o',ms=2,mec=c,c='w',mew=2)\n",
    "\n",
    "\n",
    "                pl.plot([mi+(0.25*hi),mi+(0.25*hi)],[l,u],c='w',lw=2)\n",
    "\n",
    "                pl.plot([mi+(0.25*hi),mi+(0.25*hi)],[l,u],c=c,lw=1)\n",
    "\n",
    "\n",
    "                pl.xlim(-0.5,len(this_roi_comb_order)-0.5)  \n",
    "                pl.ylim(-0.5,1)\n",
    "                pl.yticks(np.linspace(-0.5,1,4))\n",
    "    #         pl.xticks([0,stim_radius])\n",
    "\n",
    "        sn.despine(offset=2)\n",
    "\n",
    "        pl.xticks(range(len(this_roi_comb_order)),this_roi_comb_order,rotation=45)\n",
    "\n",
    "        pl.tight_layout()\n",
    "        f.savefig(os.path.join(fig_dir,'%s_consistency_%s_corrs.pdf'%(name,param)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hemis = ['']#,'left_','right_']\n",
    "for comp,name in zip([['183_1','183_2'],['181','182']],['across_runs','across_subjects']):\n",
    "    pl.close('all')\n",
    "\n",
    "    for param in ['ecc','ang','rfsize']:#,'x','y']:\n",
    "        \n",
    "        f = pl.figure(figsize=(1,1))\n",
    "        s = f.add_subplot(111)\n",
    "        pl.axhline(0,color='k')\n",
    "\n",
    "\n",
    "        this_roi_comb_order = ['OMV','VIIb','VIIIb']\n",
    "\n",
    "\n",
    "        for mi,mask_overall in enumerate(this_roi_comb_order):\n",
    "\n",
    "            for hi,hemi in enumerate(hemis):\n",
    "                \n",
    "                mask = hemi + mask_overall\n",
    "\n",
    "                # determine this roi mask\n",
    "                roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "                for subroi in roi_combs[mask]:\n",
    "                    roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "                # determine roi color\n",
    "                c = get_roi_color(mask)  \n",
    "\n",
    "                rfsize1 = data[comp[0]][roimask,dims['rfsize']]\n",
    "                ecc1 = data[comp[0]][roimask,dims['ecc']]\n",
    "                ang1 = data[comp[0]][roimask,dims['ang']]\n",
    "                xs1 = np.cos(np.radians(np.ravel(ang1))) * np.ravel(ecc1)\n",
    "                ys1 = np.sin(np.radians(np.ravel(ang1))) * np.ravel(ecc1)\n",
    "\n",
    "                rfsize2 = data[comp[1]][roimask,dims['rfsize']]\n",
    "                ecc2= data[comp[1]][roimask,dims['ecc']]\n",
    "                ang2 = data[comp[1]][roimask,dims['ang']]\n",
    "\n",
    "                xs2 = np.cos(np.radians(np.ravel(ang2))) * np.ravel(ecc2)\n",
    "                ys2 = np.sin(np.radians(np.ravel(ang2))) * np.ravel(ecc2)\n",
    "\n",
    "\n",
    "                if param == 'x':\n",
    "                    data1 = xs1\n",
    "                    data2 = xs2\n",
    "                    plotrange = [-stim_radius,stim_radius]\n",
    "                elif param == 'y':\n",
    "                    data1 = ys1\n",
    "                    data2 = ys2\n",
    "                    plotrange = [-stim_radius,stim_radius]\n",
    "\n",
    "                elif param == 'ecc':\n",
    "                    data1 = ecc1\n",
    "                    data2 = ecc2\n",
    "                    plotrange = [0,stim_radius]\n",
    "\n",
    "                elif param == 'ang':\n",
    "                    data1 = ang1\n",
    "                    data2 = ang2\n",
    "                    plotrange = [0,360]\n",
    "\n",
    "                elif param == 'rfsize':\n",
    "                    data1 = rfsize1\n",
    "                    data2 = rfsize2\n",
    "                    plotrange = [0,20]\n",
    "\n",
    "\n",
    "                v = (~np.isnan(data1))*(~np.isnan(data2))\n",
    "\n",
    "                data1 = data1[v]\n",
    "                data2 = data2[v]\n",
    "\n",
    "                if param == 'ang':\n",
    "    #                 r,p,l,u = pearsonr_ci(data1,data2)\n",
    "                    r,p,l,u = bootstrap_correlation(np.radians(data1),np.radians(data2),circ=True)\n",
    "                    print(mask,r,p)\n",
    "                else:\n",
    "                    r,p,l,u = bootstrap_correlation(data1,data2)\n",
    "    #                 r,(l,u)= pc.corrcc(np.radians(data1),np.radians(data2),ci=0.95,bootstrap_iter=int(1e4))\n",
    "\n",
    "                if np.arctanh(l) > 0:\n",
    "                    pl.plot(mi+(0.25*hi),r,'o',ms=4/len(hemis),mec=c,c=c,mew=2)\n",
    "                else:\n",
    "                    pl.plot(mi+(0.25*hi),r,'o',ms=4/len(hemis),mec=c,c='w',mew=2)\n",
    "\n",
    "\n",
    "                pl.plot([mi+(0.25*hi),mi+(0.25*hi)],[l,u],c='w',lw=2)\n",
    "\n",
    "                pl.plot([mi+(0.25*hi),mi+(0.25*hi)],[l,u],c=c,lw=1)\n",
    "\n",
    "\n",
    "                pl.xlim(-0.5,len(this_roi_comb_order)-0.5)  \n",
    "                pl.ylim(-0.5,1)\n",
    "                pl.yticks(np.linspace(-0.5,1,4))\n",
    "\n",
    "        sn.despine(offset=2)\n",
    "\n",
    "        pl.xticks(range(len(this_roi_comb_order)),this_roi_comb_order,rotation=45)\n",
    "\n",
    "        pl.tight_layout()\n",
    "        f.savefig(os.path.join(fig_dir,'%s_consistency_%s_corrs.pdf'%(name,param)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# tSNR - R2 correlations\n",
    "\n",
    "We want to find out if the amount of explained variance in the pRF model is correlated to the temporal signal to noise ratio (tSNR). The tSNR is computed in notebook 5_tSNR by deviding the mean of the signal by the standard deviation. I created these tSNRS for every voxel for every run type (bars, wedges, rings), and than averaged the resulting tSNRs over runs. \n",
    "\n",
    "In order to match the unthresholded R2s to the tSNRs, we need the original data, which is saved by subject id (and not subject rank). Therefore, we first need to uncover these ranks again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def create_roi_mask():\n",
    "    \"\"\"\n",
    "    This mask deselects all voxels outside the retinotopic\n",
    "    clusters identified in the average subject.\n",
    "    \"\"\"\n",
    "    \n",
    "    roi_fn = os.path.join(resource_dir,'volume_masks','cerebellum_retmaps.nii')\n",
    "    maskimg = nb.load(roi_fn)\n",
    "    roimask = (maskimg.get_data()>0) # valid voxels  \n",
    "\n",
    "    return roimask\n",
    "\n",
    "def determine_best_subjects():\n",
    "    \n",
    "    \"\"\"\n",
    "    This function determines the best subjects based on the median\n",
    "    r-squared within the retinotopic clusters defined in the average subject.\n",
    "    \"\"\"\n",
    "\n",
    "    roimask = create_roi_mask()\n",
    "\n",
    "    all_r2 = []\n",
    "    for sj in range(181):      \n",
    "\n",
    "        # load the prf results nifti\n",
    "        fn = os.path.join(data_dir,'all_subjects','prfresults_subject_%d.dscalar_data_sub.nii.gz'%sj)\n",
    "        img = nb.load(fn)\n",
    "        data = img.get_data()\n",
    "\n",
    "        # mask r2s with roi mask:\n",
    "        r2s = np.ravel(data[roimask,dims['r2']])    \n",
    "        # and get median within these voxels\n",
    "        all_r2.append(np.nanmedian(r2s))\n",
    "\n",
    "    # sort subjects based on median r2\n",
    "    best_subjects = np.argsort(all_r2)[::-1]\n",
    "\n",
    "    return best_subjects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "subids = determine_best_subjects()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# eccen-size plot\n",
    "this_roi_comb_order = ['all_rois','all_voxels','OMV','VIIb','VIIIb']#VIIIb','lVIIIb','VIIb','mOMV','lOMV']\n",
    "\n",
    "f = pl.figure(figsize=(2*len(this_roi_comb_order),2))\n",
    "\n",
    "for mi,mask in enumerate(this_roi_comb_order):\n",
    "\n",
    "    s = f.add_subplot(1,len(this_roi_comb_order),mi+1)\n",
    "    pl.title(mask)\n",
    "    # determine this roi mask\n",
    "    if mask == 'all_voxels':\n",
    "        roimask = cer_mask.astype(bool)\n",
    "    else:\n",
    "        roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "        for subroi in roi_combs[mask]:\n",
    "            roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "    # determine ecc limit and roi color\n",
    "    c = get_roi_color(mask)  \n",
    "    \n",
    "    all_tSNRs = []\n",
    "    all_r2s = []\n",
    "    for sj in tqdm(range(181)):\n",
    "        \n",
    "        # conver subject rank (sj) back to subject id (si) here:\n",
    "        si = subids[sj]\n",
    "        this_sj_data = nb.load(os.path.join(data_dir,'all_subjects','prfresults_subject_%d.dscalar_data_sub.nii.gz'%si)).get_data()\n",
    "        # tsnrs are saved by subject rank, so we can use the rank (sj) here:\n",
    "        tsnrs = nb.load(os.path.join(data_dir,'timeseries','sub_rank_%d_tsnr_avg_over_tasks.nii'%sj)).get_data()[roimask]\n",
    "        r2s = this_sj_data[roimask,dims['r2']]\n",
    "\n",
    "        all_r2s.append(np.nanmedian(r2s))\n",
    "        all_tSNRs.append(np.nanmedian(tsnrs))\n",
    "        \n",
    "    rho,rho_p = sp.stats.spearmanr(np.array(all_tSNRs),np.array(all_r2s))\n",
    "    r,r_p = sp.stats.pearsonr(np.array(all_tSNRs),np.array(all_r2s))\n",
    "\n",
    "    sn.regplot(np.array(all_tSNRs),np.array(all_r2s),color=c,label='r: %.2f, p: %.2f\\nrho: %.2f, p: %.2f'%(r,r_p,rho,rho_p))\n",
    "    if mi == 0:\n",
    "        pl.xlabel('tSNR')\n",
    "        pl.ylabel('R2')\n",
    "    pl.xticks([np.min(all_tSNRs),np.max(all_tSNRs)])\n",
    "    pl.yticks([np.min(all_r2s),np.max(all_r2s)])\n",
    "    \n",
    "    pl.legend(loc='best')\n",
    "    \n",
    "    sn.despine(offset=2)\n",
    "\n",
    "pl.tight_layout()\n",
    "f.savefig(os.path.join(fig_dir,'r2_tsnr_corr.pdf'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Timeseries plot\n",
    "\n",
    "In this plot, we'll get the design matrix from the osf website, and use the pRF model that was fit on all data (wedges / rings / bars) to create a predicted fit for the average of the two bar runs (downloaded in notebook 1_export_ciftis)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Popeye imports\n",
    "import popeye.utilities as utils\n",
    "from popeye.visual_stimulus import VisualStimulus\n",
    "import popeye.css as css\n",
    "import scipy as sp\n",
    "import urllib\n",
    "import zipfile\n",
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ap_dir = os.path.join(data_dir,'apertures.zip')\n",
    "ap_url = 'https://osf.io/5sj3m/download'\n",
    "up_ap_dir = os.path.join(data_dir,'apertures')\n",
    "mat_dir = os.path.join(data_dir,'RETBARsmall.mat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### download the design matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.isfile(ap_dir):\n",
    "    print('Downloading retinotopy design matrices...')\n",
    "    urllib.urlretrieve(ap_url, ap_dir)  \n",
    "    print('downloading done!')\n",
    "else:\n",
    "    print('dms already downloaded')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### unpacking the dms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.isdir(up_ap_dir):\n",
    "    print('unpacking dms')\n",
    "    zip_ref = zipfile.ZipFile(ap_dir, 'r')\n",
    "    zip_ref.extractall(ap_dir[:-4])\n",
    "    zip_ref.close()\n",
    "    # remove zip file\n",
    "    os.remove(ap_dir)\n",
    "else:\n",
    "    print('dms already unpacked')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### move and remove files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.isfile(mat_dir):\n",
    "    os.rename(os.path.join(up_ap_dir,'apertures','RETBARsmall.mat'),mat_dir)\n",
    "    shutil.rmtree(up_ap_dir)\n",
    "else:\n",
    "    print('dm ret file already moved')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### now load in file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with h5py.File(mat_dir, 'r') as mat:\n",
    "    visual_dm = mat['stim'].value.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### convert file to popey visual stimulus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create stimulus design and define model\n",
    "# ---------------------------------------\n",
    "stimulus = VisualStimulus(  stim_arr = visual_dm,\n",
    "                            viewing_distance = 102, \n",
    "                            screen_width = 29,\n",
    "                            scale_factor = 1,\n",
    "                            tr_length = 1.0,\n",
    "                            dtype = np.short)\n",
    "\n",
    "model_func = css.CompressiveSpatialSummationModel(  stimulus = stimulus,\n",
    "                                                    hrf_model = utils.spm_hrf)\n",
    "model_func.hrf_delay = 0\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### load in the averaged timecourse data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Get timeseries data\n",
    "# ---------------------------------------\n",
    "timeseries_fn = os.path.join(data_dir,'timeseries','avg_prf_timeseries.nii.gz')\n",
    "timeseries_data = nb.load(timeseries_fn).get_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# zscore timeseries\n",
    "avg = np.nanmean(timeseries_data,axis=-1)[:,:,:,np.newaxis]\n",
    "std = np.nanstd(timeseries_data,axis=-1)[:,:,:,np.newaxis]\n",
    "timeseries_data = (timeseries_data-avg)/std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### setup function that creates prediction and rescales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def create_prediction(voxel_idx,roimask):\n",
    "\n",
    "    # get prf parameters\n",
    "    ang = np.radians(avgdata[roimask,dims['ang']][voxel_idx])\n",
    "    ecc = avgdata[roimask,dims['ecc']][voxel_idx]\n",
    "    size = avgdata[roimask,dims['rfsize']][voxel_idx]\n",
    "\n",
    "    # get ts data\n",
    "    these_ts = timeseries_data[roimask]\n",
    "    \n",
    "    # convert to cartesian\n",
    "    x = ecc * np.cos(ang)\n",
    "    y = ecc * np.sin(ang)\n",
    "\n",
    "    # this is the n used in the paper\n",
    "    n = 0.05\n",
    "\n",
    "    # calculate prediction\n",
    "    pred = model_func.generate_prediction(x,y,size,n,1,0)\n",
    "    \n",
    "    # refit baseline and amp parameters (since fitted on different data)\n",
    "    dm = np.mat(np.vstack([np.ones_like(pred),pred])).T\n",
    "    t = these_ts[voxel_idx]\n",
    "    intercept,slope = np.array(np.linalg.pinv(dm.T * dm) * dm.T * np.mat(t[:,np.newaxis]))     \n",
    "\n",
    "    # scale prediction\n",
    "    p =(intercept[0]+pred*slope[0])\n",
    "\n",
    "    # compute new r2\n",
    "    r2 = (sp.stats.pearsonr(p,t)[0]**2)*100\n",
    "    \n",
    "    return p,r2\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine best voxels for these data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "these_r2s = {}\n",
    "ecc_bands = [[0,2],[2,4],[4,6],[6,8]]\n",
    "for ecc_band in ecc_bands:\n",
    "    these_r2s[str(ecc_band)] = {}\n",
    "    \n",
    "    # select data from this ecc band\n",
    "    eccs = avgdata[:,:,:,dims['ecc']]\n",
    "    valid_eccs = (eccs>ecc_band[0])*(eccs<ecc_band[1])\n",
    "    \n",
    "    for mi,mask in enumerate(roi_comb_order):\n",
    "        print('creating predictions for ecc_band %s, mask %s'%(ecc_band,mask))\n",
    "        # determine this roi mask\n",
    "        roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "        for subroi in roi_combs[mask]:\n",
    "            roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "        roimask *= valid_eccs\n",
    "\n",
    "        # apply mask to timeseries\n",
    "        these_ts = timeseries_data[roimask]\n",
    "    \n",
    "        # now create a prediction for every voxel\n",
    "        # on the averaged data and save the r2s\n",
    "        r2s = []\n",
    "        for v in range(roimask.sum()):\n",
    "\n",
    "            p,r2 = create_prediction(v,roimask)\n",
    "\n",
    "            # recompute r2 on this prediction and data\n",
    "            r2s.append(r2)\n",
    "            \n",
    "        these_r2s[str(ecc_band)][mask] = r2s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### now we can create plot for the best voxels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create plot\n",
    "# ---------------------------------------\n",
    "n_vox = 5\n",
    "for ecc_band in ecc_bands:\n",
    "    \n",
    "    # select data from this ecc band\n",
    "    eccs = avgdata[:,:,:,dims['ecc']]\n",
    "    valid_eccs = (eccs>ecc_band[0])*(eccs<ecc_band[1])   \n",
    "    \n",
    "    for mi,mask in enumerate(roi_comb_order):\n",
    "\n",
    "        # determine this roi mask\n",
    "        roimask = np.zeros_like(cer_retmaps).astype(bool)\n",
    "        for subroi in roi_combs[mask]:\n",
    "            roimask[(cer_retmaps==subroi)] = True\n",
    "\n",
    "        roimask *= valid_eccs\n",
    "\n",
    "        # apply mask to timeseries\n",
    "        these_ts = timeseries_data[roimask]\n",
    "\n",
    "        # now determine best voxel based on fits to this data\n",
    "        r2s = these_r2s[str(ecc_band)][mask]\n",
    "        best_voxels = np.argsort(r2s)[::-1][:n_vox]\n",
    "\n",
    "        # and create plots for these voxels\n",
    "        for best_voxel in best_voxels:\n",
    "\n",
    "            # get this timeseries\n",
    "            timeseries = these_ts[best_voxel]\n",
    "\n",
    "            pred,r2 = create_prediction(best_voxel,roimask)\n",
    "            \n",
    "            # generate pRF\n",
    "            res = 501\n",
    "\n",
    "            # get pRF parameters for visualization\n",
    "            ang = np.radians(avgdata[roimask,dims['ang']][best_voxel])\n",
    "            ecc = avgdata[roimask,dims['ecc']][best_voxel]\n",
    "            size = avgdata[roimask,dims['rfsize']][best_voxel]\n",
    "            \n",
    "            x = ecc * np.cos(ang)\n",
    "            y = ecc * np.sin(ang)\n",
    "\n",
    "            f=pl.figure(figsize=(4,1.25))            \n",
    "            \n",
    "            # first plot the prf parameters\n",
    "            s = f.add_subplot(141)\n",
    "            pl.text(0,0,'ecc: %.2f\\nangle: %.2f\\nsize: %.2f\\nR2: %.2f'%(ecc,ang,size,r2),\n",
    "                    horizontalalignment='right',verticalalignment='center')            \n",
    "            pl.xlim(-5,5)\n",
    "            pl.ylim(-5,5)\n",
    "            pl.axis('off')\n",
    "            \n",
    "        \n",
    "            # now create a pRF visualization plot\n",
    "            s = f.add_subplot(142,aspect='equal')\n",
    "\n",
    "            c = get_roi_color(mask)\n",
    "\n",
    "            # crosshair\n",
    "            pl.axhline(0,lw=0.5,color='k')\n",
    "            pl.axvline(0,lw=0.5,color='k')\n",
    "            \n",
    "            # for prf size\n",
    "            s.add_artist(pl.Circle((x,y),size, color='r',fill=False,alpha=0.25))\n",
    "            # prf center\n",
    "            pl.plot(x,y,'o',color='r',ms=3,mec='w',mew=1,alpha=1)\n",
    "\n",
    "            # plot properties\n",
    "            pl.xlim(-stim_radius,stim_radius)\n",
    "            pl.ylim(-stim_radius,stim_radius) \n",
    "            pl.xticks([-stim_radius,0,stim_radius])\n",
    "            pl.yticks([-stim_radius,0,stim_radius])\n",
    "            pl.xlabel('visual field x (dva)')\n",
    "            pl.ylabel('visual field y (dva)')\n",
    "            sn.despine(offset=2)\n",
    "            \n",
    "            # and finally the pRF timeseries prediction\n",
    "            s = f.add_subplot(122)\n",
    "    \n",
    "            # timeseries\n",
    "            pl.plot(timeseries,'o--',color='k',ms=1.5)\n",
    "            # prediction\n",
    "            pl.plot(pred,color='r',lw=1.5)\n",
    "\n",
    "            # this is the gray/white shading to indicate stimulus design\n",
    "            ylims=s.get_ylim()\n",
    "            pl.fill_between([0,16],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='w')\n",
    "            pl.fill_between([16,48],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='k',alpha=0.1)\n",
    "            pl.fill_between([48,80],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='w')\n",
    "            pl.fill_between([80,112],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='k',alpha=0.1)\n",
    "            pl.fill_between([112,144],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='w')\n",
    "            pl.fill_between([144,156],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='k',alpha=0.1)\n",
    "            pl.fill_between([156,188],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='w')\n",
    "            pl.fill_between([188,220],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='k',alpha=0.1)\n",
    "            pl.fill_between([220,252],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='w')\n",
    "            pl.fill_between([252,284],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='k',alpha=0.1) \n",
    "            pl.fill_between([284,300],[ylims[0],ylims[0]],[ylims[1],ylims[1]],color='w')             \n",
    "                \n",
    "            # plot properties\n",
    "            sn.despine(offset=2)\n",
    "            pl.xlabel('time (s)')\n",
    "            pl.ylabel('BOLD (z-score)')\n",
    "\n",
    "            # save figure\n",
    "            pl.tight_layout()\n",
    "            f.savefig(os.path.join(fig_dir,'pred_%s_eccband_%s_v_%d.pdf'%(mask,ecc_band,best_voxel)))\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [analysis]",
   "language": "python",
   "name": "Python [analysis]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}