create_figures.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import pickle\n",
    "import nibabel\n",
    "import cortex\n",
    "from IPython.core.debugger import set_trace\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import warnings\n",
    "import subprocess\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "if not os.path.exists(\"figures\"):\n",
    "    os.mkdir(\"figures\")\n",
    "    \n",
    "if not os.path.exists(\"figures/r2plus_figures\"):\n",
    "    os.mkdir(\"figures/r2plus_figures\")\n",
    "\n",
    "if not os.path.exists(\"figures/r2plus_diff_figures\"):\n",
    "    os.mkdir(\"figures/r2plus_diff_figures\")\n",
    "    \n",
    "if not os.path.exists(\"figures/sig_figures\"):\n",
    "    os.mkdir(\"figures/sig_figures\")\n",
    "    \n",
    "if not os.path.exists(\"figures/roi_figures\"):\n",
    "    os.mkdir(\"figures/roi_figures\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_pattern = \"{base}_{view}_{surface}.png\"\n",
    "_combine = lambda a,b: ( lambda c: [c, c.update(b)][0] )(dict(a))\n",
    "_tolists = lambda p: {k:[v] for k,v in p.items()}\n",
    "\n",
    "def save_3d_views(data, root, base_name, list_views =['lateral'],list_surfaces = ['fiducial'], with_labels = False,\n",
    "                  size=(1024*4, 768*4), trim=True):\n",
    "    \"\"\"Saves 3D views of `data` in and around `root` under multiple specifications. Needs to be run\n",
    "       on a system with a display (will launch webgl viewer)\n",
    "    data: a pycortex volume\n",
    "    root: directory where things should be saved\n",
    "    base_name: base name for images\n",
    "    list_views: which views do you want? choices are: lateral, lateral_left, lateral_right,\n",
    "               medial, front, back,top, bottom\n",
    "    list_surfaces: what surfaces do you want? choices are inflated, flatmap, fiducial\n",
    "    with_labels: show ROI labels?\n",
    "    size: size of produced image (before trimming)\n",
    "    trim: whether to trim\n",
    "    returns filenames: a dict of the produced image paths\n",
    "    \"\"\"\n",
    "    # Create root dir?\n",
    "    if not os.path.exists(root):\n",
    "        os.mkdir(root)\n",
    "\n",
    "    # Create viewer\n",
    "    if with_labels:\n",
    "        labels_visible=('rois',)\n",
    "    else:\n",
    "        labels_visible=()\n",
    "    handle = cortex.webgl.show(data,labels_visible=labels_visible)\n",
    "\n",
    "    time.sleep(5.0)\n",
    "\n",
    "\n",
    "    basic = dict()#radius=400)#projection=['orthographic'], #radius=260, visL=True, visR=True)\n",
    "\n",
    "    views = dict(lateral=dict(altitude=80.5, azimuth=181, pivot=180.5, radius = 700),#, specularity = 0),\n",
    "                 lateral_left=dict(altitude=90.5, azimuth=90.5, pivot=0.5),\n",
    "                 lateral_right=dict(altitude=90.5, azimuth=270.5, pivot=0.5),\n",
    "                 medial=dict(altitude=90.5, azimuth=0.5, pivot=180.5),\n",
    "                 front=dict(altitude=90.5, azimuth=0, pivot=0),\n",
    "                 back=dict(altitude=90.5, azimuth=181, pivot=0),\n",
    "                 top=dict(altitude=0, azimuth=180, pivot=0),\n",
    "                 bottom=dict(altitude=180, azimuth=0, pivot=0)\n",
    "                )\n",
    "\n",
    "    surfaces = dict(inflated=dict(unfold= 0.5) ,\n",
    "                    flatmap=dict(unfold= 1) ,\n",
    "                    fiducial=dict(unfold= 0)\n",
    "                   )\n",
    "\n",
    "    param_dict = dict(unfold = 'surface.{subject}.unfold',\n",
    "                      altitude = 'camera.altitude',\n",
    "                      azimuth = 'camera.azimuth',\n",
    "                      pivot = 'surface.{subject}.pivot',\n",
    "                      radius = 'camera.radius')\n",
    "#                      specularity = 'surface.{subject}.specularity') # unknown parameter\n",
    "\n",
    "\n",
    "    # Save views!\n",
    "    filenames = dict([(key, dict()) for key in surfaces.keys()])\n",
    "\n",
    "    for view in list_views:\n",
    "        # copy proper parameters with new names\n",
    "        vparams = dict([(param_dict[k], v) for k,v in views[view].items()])\n",
    "        for surf in list_surfaces:\n",
    "            # copy proper parameters with new names\n",
    "            sparams = dict([(param_dict[k], v) for k,v in surfaces[surf].items()])\n",
    "            # Combine basic, view, and surface parameters\n",
    "            params = _combine(_combine(basic, vparams), sparams)\n",
    "\n",
    "            # Set the view\n",
    "            handle._set_view(**_tolists(params))\n",
    "            time.sleep(5.5)\n",
    "\n",
    "            # Save image, store filename\n",
    "            filename = file_pattern.format(base=base_name, view=view, surface=surf)\n",
    "            filenames[surf][view] = filename\n",
    "            # filenames.append(filename)\n",
    "\n",
    "            output_path = os.path.join(root, filename)\n",
    "            print(output_path)\n",
    "            handle.getImage(output_path, size)\n",
    "            \n",
    "            while not os.path.exists(output_path):\n",
    "                pass\n",
    "            time.sleep(0.1)\n",
    "\n",
    "            # Trim edges?\n",
    "            if trim:\n",
    "                # Wait for browser to dump file\n",
    "                while not os.path.exists(output_path):\n",
    "                    pass\n",
    "\n",
    "                time.sleep(0.5)\n",
    "                subprocess.call([\"convert\", \"-trim\", output_path, output_path])\n",
    "\n",
    "    # Close the window!\n",
    "    try:\n",
    "        handle.close()\n",
    "    except:\n",
    "        print('Could not close viewer')\n",
    "\n",
    "    return filenames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sqrt R2 figures\n",
    "def sqrt_r2(vec):\n",
    "    vec2=np.zeros_like(vec)\n",
    "    vec2[vec>0] = vec[vec>0]\n",
    "    vec2 = np.sqrt(vec2)\n",
    "    return vec2\n",
    "\n",
    "predictions_dir = \"predictions_mni/\"\n",
    "\n",
    "mask = pickle.load(open(\"mni_mask.pkl\",\"rb\"))\n",
    "vols = {}\n",
    "all_subjects = [\"F\",\"G\",\"H\",\"I\",\"J\",\"K\",\"L\",\"M\",\"N\"]\n",
    "\n",
    "all_features = [\"punct_final\",\n",
    "                \"all_complexity_metrics_punct\",\n",
    "                \"pos_dep_tags_all_complexity_metrics\", \n",
    "                \"aggregated_contrege_comp_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_incontrege_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_bert_PCA_dims_15_contrege_incomp_pos_dep_tags_all_complexity_metrics\"]\n",
    "\n",
    "for feat in all_features:\n",
    "    x = 0\n",
    "    for sub in all_subjects:\n",
    "        p = np.load(predictions_dir + \"{}/{}_r2s.npy\".format(feat, sub))\n",
    "        if x is 0:\n",
    "            x = np.zeros(p.shape)\n",
    "        x += p/len(all_subjects)\n",
    "    vols[feat] = x\n",
    "\n",
    "v = {}\n",
    "for feat in vols:\n",
    "    v[feat] = cortex.Volume(sqrt_r2(vols[feat]), 'MNI','atlas',mask=mask, vmin = 0, vmax = 0.15, cmap=\"Reds\")\n",
    "    save_3d_views(v[feat], 'figures/r2plus_figures/',\"r2plus_{}\".format(feat), trim=True, list_surfaces = ['inflated'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sqrt R2 figures\n",
    "def sqrt_r2(vec):\n",
    "    vec2=np.zeros_like(vec)\n",
    "    vec2[vec>0] = vec[vec>0]\n",
    "    vec2 = np.sqrt(vec2)\n",
    "    return vec2\n",
    "\n",
    "predictions_dir = \"predictions_mni/\"\n",
    "\n",
    "mask = pickle.load(open(\"mni_mask.pkl\",\"rb\"))\n",
    "vols = {}\n",
    "all_subjects = [\"F\",\"G\",\"H\",\"I\",\"J\",\"K\",\"L\",\"M\",\"N\"]\n",
    "\n",
    "all_feature_pairs = [\n",
    "                     (\"all_complexity_metrics_punct\", \"punct_final\"), # (CM + PU) - (PU)\n",
    "                     (\"pos_dep_tags_all_complexity_metrics\", \"all_complexity_metrics_punct\"), # since PD already contains PU, this tests (PD + CM + PU) - (CM + PU)\n",
    "                     (\"aggregated_contrege_comp_pos_dep_tags_all_complexity_metrics\", \"pos_dep_tags_all_complexity_metrics\"), # (CC + PD + CM + PU) - (PD + CM + PU)\n",
    "                     (\"aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics\", \"pos_dep_tags_all_complexity_metrics\"), # (CI + PD + CM + PU) - (PD + CM + PU)\n",
    "                     (\"aggregated_incontrege_pos_dep_tags_all_complexity_metrics\", \"pos_dep_tags_all_complexity_metrics\"), # (INC + PD + CM + PU) - (PD + CM + PU)\n",
    "                     (\"aggregated_bert_PCA_dims_15_contrege_incomp_pos_dep_tags_all_complexity_metrics\", \"aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics\") # (BERT + CI + PD + CM + PU) - (CI + PD + CM + PU)\n",
    "                    ]\n",
    "\n",
    "all_features = [\"punct_final\", \n",
    "                \"all_complexity_metrics_punct\",\n",
    "                \"pos_dep_tags_all_complexity_metrics\", \n",
    "                \"aggregated_contrege_comp_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_incontrege_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_bert_PCA_dims_15_contrege_incomp_pos_dep_tags_all_complexity_metrics\"]\n",
    "\n",
    "all_base_vols = {}\n",
    "for feat in all_features:\n",
    "    x = 0\n",
    "    for sub in all_subjects:\n",
    "        p = np.load(predictions_dir + \"{}/{}_r2s.npy\".format(feat, sub))\n",
    "        if x is 0:\n",
    "            x = np.zeros(p.shape)\n",
    "        x += p/len(all_subjects)\n",
    "    all_base_vols[feat] = x\n",
    "\n",
    "v = {}\n",
    "for pair in all_feature_pairs:\n",
    "    diff = sqrt_r2(all_base_vols[pair[0]]) - sqrt_r2(all_base_vols[pair[1]])\n",
    "    v[\"{}_diff_{}\".format(pair[0], pair[1])] = cortex.Volume(diff, 'MNI','atlas',mask=mask, vmin = 0, vmax = 0.05, cmap=\"Reds\")\n",
    "    save_3d_views(v[\"{}_diff_{}\".format(pair[0], pair[1])], 'figures/r2plus_diff_figures/',\"r2plus_diff_{}\".format(\"{}_diff_{}\".format(pair[0], pair[1])), trim=True, list_surfaces = ['inflated'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sig corrected figures\n",
    "predictions_dir = \"predictions_mni/\"\n",
    "\n",
    "mask = pickle.load(open(\"mni_mask.pkl\",\"rb\"))\n",
    "vols = {}\n",
    "all_subjects = [\"F\",\"G\",\"H\",\"I\",\"J\",\"K\",\"L\",\"M\",\"N\"]\n",
    "\n",
    "all_features = [\"punct_final\",\n",
    "                \"all_complexity_metrics_punct_diff_punct_final\",\n",
    "                \"pos_dep_tags_all_complexity_metrics_diff_all_complexity_metrics_punct\",\n",
    "                \"aggregated_contrege_comp_pos_dep_tags_all_complexity_metrics_diff_pos_dep_tags_all_complexity_metrics\", \n",
    "                \"aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics_diff_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_incontrege_pos_dep_tags_all_complexity_metrics_diff_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_bert_PCA_dims_15_contrege_incomp_pos_dep_tags_all_complexity_metrics_diff_aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics\"]\n",
    "\n",
    "for feat in all_features:\n",
    "    x = 0\n",
    "    for sub in all_subjects:\n",
    "        p = 0\n",
    "        if os.path.exists(predictions_dir + \"{}/{}_sig_bootstrap_group_corrected.npy\".format(feat, sub)):\n",
    "            p = np.load(predictions_dir + \"{}/{}_sig_bootstrap_group_corrected.npy\".format(feat, sub))\n",
    "        elif os.path.exists(predictions_dir + \"{}/{}_sig_group_corrected.npy\".format(feat, sub)):\n",
    "            p = np.load(predictions_dir + \"{}/{}_sig_group_corrected.npy\".format(feat, sub))\n",
    "        if x is 0:\n",
    "            x = np.zeros(p.shape)\n",
    "        x[p!=0] += 1\n",
    "    vols[feat] = x\n",
    "\n",
    "v = {}\n",
    "for feat in vols:\n",
    "    v[feat] = cortex.Volume(vols[feat], 'MNI','atlas',mask=mask, vmin = 0, vmax = 5, cmap = \"J4\")\n",
    "    save_3d_views(v[feat], 'figures/sig_figures/',\"sig_{}\".format(feat), trim=True, list_surfaces = ['inflated'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# blank plot of all ROIs\n",
    "roi_file_name = os.path.join('LangParcels_n220_LH.hdr')\n",
    "roi_dat = nibabel.load(roi_file_name)\n",
    "roi_mat = roi_dat.get_data().T\n",
    "mask_narrow_mni = np.zeros((91,109,91))\n",
    "mask_narrow_mni[11:-11,7:-7,6:-6] = roi_mat\n",
    "# flip\n",
    "mask_narrow_mni[:,:,:45] = mask_narrow_mni[:,:,46:][:,:,::-1]\n",
    "plt.hist(mask_narrow_mni[mask_narrow_mni>0],100)\n",
    "print(np.unique(mask_narrow_mni))\n",
    "for i in range(1,7):\n",
    "    mask_narrow_mni[ (mask_narrow_mni>i-0.5)*(mask_narrow_mni<i+0.5)] = i\n",
    "print(np.unique(mask_narrow_mni))\n",
    "plt.hist(mask_narrow_mni[mask_narrow_mni>0],100)\n",
    "plt.figure()\n",
    "np.save(\"processed_fedorenko_masks_MNI\", mask_narrow_mni)\n",
    "vols_parcels_narrow = cortex.Volume(mask_narrow_mni,'MNI','atlas_2mm', vmin = 0, vmax = 6, cmap=\"gist_ncar_r\")\n",
    "save_3d_views(vols_parcels_narrow, 'figures/',\"blank_plot_of_rois\", trim=True, list_surfaces = ['inflated'], with_labels=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ROI analysis\n",
    "np.random.seed(97)\n",
    "predictions_dir = 'predictions/'\n",
    "all_subjects = [\"F\",\"G\",\"H\",\"I\",\"J\",\"K\",\"L\",\"M\",\"N\"]\n",
    "\n",
    "mask_sub = np.load('mask_sub_fedorenko_HP_narrow_mirrored.npy', allow_pickle=True, fix_imports=True)[()] # because it was saved with np.save\n",
    "\n",
    "mask_num = dict([(s,mask_sub[s][np.load(\"masks/{}.npy\".format(s))])\n",
    "                 for s in all_subjects])\n",
    "\n",
    "roi_names = ['PTL',\n",
    "             'ATL',\n",
    "             'AG',\n",
    "             'IFG',\n",
    "             'MFG',\n",
    "             'IFGorb']\n",
    "\n",
    "all_features = [\"all_complexity_metrics_punct_diff_punct_final\",\n",
    "                \"pos_dep_tags_all_complexity_metrics_diff_all_complexity_metrics_punct\",\n",
    "                \"aggregated_contrege_comp_pos_dep_tags_all_complexity_metrics_diff_pos_dep_tags_all_complexity_metrics\", \n",
    "                \"aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics_diff_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_incontrege_pos_dep_tags_all_complexity_metrics_diff_pos_dep_tags_all_complexity_metrics\",\n",
    "                \"aggregated_bert_PCA_dims_15_contrege_incomp_pos_dep_tags_all_complexity_metrics_diff_aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics\"]\n",
    "\n",
    "metric = \"sig_bootstrap_group_corrected\"\n",
    "\n",
    "n_rois = len(roi_names)\n",
    "NFS = len(all_features)\n",
    "NS = len(all_subjects)\n",
    "percent_sig = []\n",
    "np.zeros((NFS, NS, n_rois))\n",
    "\n",
    "all_names = [\"{CM, PU} - {PU}\", \n",
    "             \"{PD, CM, PU} - {CM, PU}\", \n",
    "             \"{CC, PD, CM, PU} - {PD, CM, PU}\", \n",
    "             \"{CI, PD, CM, PU} - {PD, CM, PU}\",\n",
    "             \"{INC, PD, CM, PU} - {PD, CM, PU}\",\n",
    "             \"{BERT, CI, PD, CM, PU} - {CI, PD, CM, PU}\"]\n",
    "\n",
    "for i_s, sub in enumerate(all_subjects):\n",
    "    mask = np.load(\"masks/{}.npy\".format(sub))\n",
    "    for i, roi in enumerate(roi_names):        \n",
    "        roi_mask = np.where(mask_num[sub]==(i+1))[0]\n",
    "        \n",
    "        for ifs, feature in enumerate(all_features):\n",
    "            feature_metrics = np.load(\"{}{}/{}_{}.npy\".format(predictions_dir, feature, sub, metric))\n",
    "            if len(feature_metrics.shape) > 1 and feature_metrics.shape[1] > 1:\n",
    "                feature_metrics = feature_metrics[:,0]\n",
    "            \n",
    "            data_point = {}\n",
    "            data_point[\"Test\"] = all_names[ifs]\n",
    "            data_point[\"Feature Added\"] = all_names[ifs].split(\",\")[0][1:]\n",
    "            data_point[\"Subject\"] = sub\n",
    "            data_point[\"ROI\"] = roi\n",
    "            data_point[\"% of significant voxels\"] = np.nan_to_num(((feature_metrics[roi_mask] != 0).sum() / feature_metrics[roi_mask].shape[0]) * 100, posinf = 0, neginf = 0)\n",
    "            \n",
    "            percent_sig.append(data_point)\n",
    "            \n",
    "percent_sig = pd.DataFrame(percent_sig)\n",
    "\n",
    "fontsize = 20\n",
    "fig, ax = plt.subplots(figsize=(15,8))\n",
    "plt.rcParams.update({'font.size': fontsize})\n",
    "sns.barplot(x=\"ROI\", y=r\"% of ROI voxels with significant $R^2$ increases\", hue=\"Test\", data=percent_sig, ax=ax, ci=68)\n",
    "ax.legend(fontsize = fontsize, bbox_to_anchor=(0.5, 1.16), loc='center',frameon=False, ncol=2)\n",
    "ax.autoscale_view()\n",
    "plt.savefig(os.path.join(\"figures/roi_figures\", \"combined.png\"), bbox_inches='tight')\n",
    "plt.clf()\n",
    "\n",
    "for ROI in roi_names:\n",
    "    data_for_this_roi = percent_sig[percent_sig[\"ROI\"] == ROI]\n",
    "    fig, ax = plt.subplots(figsize=(15,8))\n",
    "    ax.set_title(ROI)\n",
    "    sns.barplot(x=\"Feature Added\", y=r\"% of ROI voxels with significant $R^2$ increases\", \\\n",
    "                data=data_for_this_roi, ci=68, ax=ax, linewidth=2.5, facecolor=(1, 1, 1, 0), errcolor=\".2\", edgecolor=\".2\", capsize=0.25)\n",
    "    sns.scatterplot(x=\"Feature Added\", y=r\"% of ROI voxels with significant $R^2$ increases\", data=data_for_this_roi, ax=ax, color=\"red\")\n",
    "    plt.savefig(os.path.join(\"figures/roi_figures\", \"{}.png\".format(ROI)), bbox_inches='tight')\n",
    "    plt.clf()"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}