loadDataProcessing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "positive-sampling",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:24.958346Z",
     "start_time": "2022-04-20T06:03:24.956007Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# load data and processing\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "medieval-fifty",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:28.181856Z",
     "start_time": "2022-04-20T06:03:25.349954Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy.linalg import norm\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy.matlib\n",
    "import pandas as pd\n",
    "\n",
    "from allensdk.api.queries.ontologies_api import OntologiesApi\n",
    "from allensdk.core.mouse_connectivity_cache import MouseConnectivityCache\n",
    "from mcmodels.core import VoxelModelCache\n",
    "\n",
    "import pickle\n",
    "\n",
    "import networkx as nx\n",
    "import plotly.offline as py\n",
    "import plotly.graph_objects as go\n",
    "import re\n",
    "\n",
    "import sklearn as sk\n",
    "# import plotly.express as px"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "breeding-pillow",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:28.877775Z",
     "start_time": "2022-04-20T06:03:28.184081Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import plotly.express as px\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42\n",
    "plt.rcParams['svg.fonttype'] = 'none'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hazardous-instrumentation",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:28.888621Z",
     "start_time": "2022-04-20T06:03:28.880108Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# load helper functions\n",
    "from helper import load_hierarchy, plot_connectivity, generate_connectivity, generate_random_connectivity\n",
    "from helper import load_celldensity, load_interneurondensity, interneuronDensityProcessing\n",
    "import modelHelper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "particular-programming",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:28.902762Z",
     "start_time": "2022-04-20T06:03:28.890479Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "filepath = './'  # change path accordingly\n",
    "libfilepath = './lib/'  # change path accordingly\n",
    "\n",
    "# load the harris hierarchy N =37; \n",
    "area_list_hier = ['ACAd', 'ACAv', 'AId', 'AIp', 'AIv', 'AUDd', 'AUDp', 'AUDpo', 'FRP', 'ILA', 'MOp', 'MOs', 'ORBl',\n",
    "         'ORBm', 'ORBvl', 'PL', 'RSPagl', 'RSPd', 'RSPv', 'SSp-bfd', 'SSp-ll', 'SSp-m', 'SSp-n',\n",
    "         'SSp-tr', 'SSp-ul', 'SSs', 'TEa', 'VISa', 'VISal', 'VISam', 'VISl', 'VISli', 'VISp',\n",
    "         'VISpl', 'VISpm', 'VISpor', 'VISrl']  # list of areas in the current model based on published paper.\n",
    "        # SSp-un and VISC are removed. \n",
    "\n",
    "hierfilename = 'mouse_harrishier_extracted_CC_published.csv'\n",
    "\n",
    "[harrisHierarchy, harrisHierarchy_df] = load_hierarchy(libfilepath, hierfilename, area_list_hier)\n",
    "harrisHierarchy_df.to_csv('results/dfHarrisHier.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "improving-lawrence",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:28.917299Z",
     "start_time": "2022-04-20T06:03:28.904692Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "harrisHierarchy_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "behavioral-edwards",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:29.580452Z",
     "start_time": "2022-04-20T06:03:29.576372Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "area_list = ['ACAd', 'ACAv', 'AId', 'AIp', 'AIv', 'AUDd', 'AUDp', 'AUDpo', 'AUDv', 'ECT', 'FRP', 'GU', 'ILA', 'MOp', 'MOs', 'ORBl',\n",
    "         'ORBm', 'ORBvl', 'PERI', 'PL', 'RSPagl', 'RSPd', 'RSPv', 'SSp-bfd', 'SSp-ll', 'SSp-m', 'SSp-n',\n",
    "         'SSp-tr', 'SSp-ul', 'SSp-un', 'SSs', 'TEa','VISC', 'VISa', 'VISal', 'VISam', 'VISl', 'VISli', 'VISp',\n",
    "         'VISpl', 'VISpm', 'VISpor', 'VISrl']\n",
    "n_areas = len(area_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "worthy-divide",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:30.418081Z",
     "start_time": "2022-04-20T06:03:30.413120Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# define divisions from Harris et al.\n",
    "div = {'somatomotor':['SSp-bfd','SSp-tr','SSp-ll','SSp-ul',\n",
    "                     'SSp-un','SSp-n','SSp-m','SSs','MOp','MOs'],\n",
    "      'visual':['VISal','VISl','VISp','VISpl',\n",
    "                'VISli','VISpor','VISrl'],\n",
    "      'medial':['VISa','VISam','VISpm','RSPagl',\n",
    "               'RSPd','RSPv'],\n",
    "      'prefrontal':['FRP','ACAd','ACAv','PL',\n",
    "                   'ILA','ORBl','ORBm','ORBvl'],\n",
    "      'lateral':['AId','AIv','AIp','VISC',\n",
    "                 'GU','TEa','PERI','ECT'],\n",
    "      'auditory':['AUDd','AUDp','AUDpo','AUDv']}\n",
    "\n",
    "div_color_list = np.array([[255, 140, 0],\n",
    "                           [169, 164, 217],\n",
    "                           [22, 37, 128],\n",
    "                           [220,20,60],\n",
    "                           [218,165,32],\n",
    "                           [91,  0,  105]])/255\n",
    "# div_color_list = {\n",
    "#     'somatomotor':[255, 140, 0],\n",
    "#       'visual':[169, 164, 217],\n",
    "#       'medial':[22, 37, 128],\n",
    "#       'prefrontal':[220,20,60],\n",
    "#       'lateral':[218,165,32],\n",
    "#       'auditory':[91,  0,  105]}\n",
    "\n",
    "div_name_list = ['somatomotor',\n",
    "                 'visual',\n",
    "                 'medial',\n",
    "                 'prefrontal',\n",
    "                 'lateral',\n",
    "                 'auditory']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "injured-alignment",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:30.882207Z",
     "start_time": "2022-04-20T06:03:30.873898Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def plotHiernotebook(hierarchy_df, dotSize, yaxis, figureSize, fontSize, saveFig, fileName):\n",
    "        # seaborn.set()\n",
    "    dfRaw = hierarchy_df.sort_values(hierarchy_df.columns[0])\n",
    "\n",
    "    fig = plt.figure(figsize = figureSize, dpi=300, facecolor='w', edgecolor='k')\n",
    "    plt.rcParams.update({'font.size': fontSize})\n",
    "    \n",
    "    xx = np.arange(len(hierarchy_df))\n",
    "    plt.plot(xx, dfRaw.values, '.-', markersize=dotSize, color='k')\n",
    "    plt.xticks(xx, dfRaw.index, rotation = 90, fontsize=fontSize*0.6)\n",
    "    plt.ylabel(yaxis) #r'$e_i/e_1$'\n",
    "    # regPV = sk.linear_model.LinearRegression()\n",
    "    # xx2 = xx[np.newaxis].T\n",
    "    # regPV.fit(xx2, dfRaw.values)\n",
    "    # print(regPV.coef_, regPV.intercept_)\n",
    "    # plt.plot(xx, xx * regPV.coef_[0][0] + regPV.intercept_[0], alpha = 0.5)\n",
    "    ax = plt.gca()\n",
    "#     [t.set_color('red') for t in ax.xaxis.get_ticklabels()]\n",
    "    for t in ax.xaxis.get_ticklabels():\n",
    "        label = t._text\n",
    "        for i in range(6):\n",
    "            if label in div[div_name_list[i]]:\n",
    "                t.set_color(div_color_list[i])\n",
    "    # plt.grid(alpha = 0.5)\n",
    "\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    if saveFig:\n",
    "        plt.savefig('figure/' + fileName, dpi=300, bbox_inches='tight',transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "romance-planning",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:31.830085Z",
     "start_time": "2022-04-20T06:03:31.272604Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plotHiernotebook(harrisHierarchy_df, 8, r'Harris Hierarchy', (5,3), 12, True, 'loadDataProcessing/harrisHierarchyPlotnotebook.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "resistant-trace",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:32.856150Z",
     "start_time": "2022-04-20T06:03:32.709261Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# load connectivity matrix\n",
    "mcc = MouseConnectivityCache(manifest_file='connectivity/mouse_connectivity_manifest.json')\n",
    "# The manifest file is a simple JSON file that keeps track of all of\n",
    "# the data that has already been downloaded onto the hard drives.\n",
    "# If you supply a relative path, it is assumed to be relative to your\n",
    "# current working directory.\n",
    "\n",
    "# grab the StructureTree instance\n",
    "structure_tree = mcc.get_structure_tree()\n",
    "\n",
    "oapi = OntologiesApi()\n",
    "\n",
    "# get the ids of all the structure sets in the tree\n",
    "# structure_set_ids = structure_tree.get_structure_sets()\n",
    "\n",
    "# download and cache the latest voxel model\n",
    "# this method returns a tuple with object types:\n",
    "# (VoxelConnectivityArray, Mask, Mask)\n",
    "cache = VoxelModelCache(manifest_file='connectivity/voxel_model_manifest.json')\n",
    "\n",
    "# extracting the cortical subnetwork. get the normalized connectvity density for ipsi side.\n",
    "normalized_connection_density = cache.get_normalized_connection_density()\n",
    "full_matrix = normalized_connection_density['ipsi']  # .loc[cortex_ids_int][cortex_ids_str]\n",
    "mat = np.array(full_matrix)\n",
    "\n",
    "# cortex\n",
    "cortex_structures = structure_tree.get_structures_by_set_id([688152357]) # first get the full info of cortex area\n",
    "# cortex_structures_refined = []\n",
    "# for s in cortex_structures:\n",
    "#     if s['id'] in full_matrix.index:\n",
    "#         cortex_structures_refined.append(s)\n",
    "cortex_ids_int = [s['id'] for s in cortex_structures]  # get the cortex area id\n",
    "cortex_ids_str = [str(s['id']) for s in cortex_structures]  \n",
    "cortex_acr = [s['acronym'] for s in cortex_structures]\n",
    "cortex_names = [s['name'] for s in cortex_structures]\n",
    "\n",
    "#thalamus\n",
    "thal_structures = structure_tree.get_structures_by_set_id([688152363])\n",
    "thal_structures_refined = [] # some of the thalmic structure not in the cconnectivity matirx\n",
    "for s in thal_structures:\n",
    "    if s['id'] in full_matrix.index:\n",
    "        thal_structures_refined.append(s)\n",
    "\n",
    "thal_ids_int = [s['id'] for s in thal_structures_refined] # store as a int\n",
    "thal_ids_str = [str(s['id']) for s in thal_structures_refined] # store as a str\n",
    "thal_acr = [s['acronym'] for s in thal_structures_refined]\n",
    "thal_names = [s['name'] for s in thal_structures_refined]\n",
    "\n",
    "\n",
    "cortical_area_number = np.size(cortex_acr)\n",
    "print(cortical_area_number)\n",
    "\n",
    "# W = mat[0:cortical_area_number, 0:cortical_area_number]  # extract the cortical areas\n",
    "\n",
    "W_ = full_matrix.loc[cortex_ids_int][cortex_ids_str] # this gives the cortico cortical connection\n",
    "W_cxth_ = full_matrix.loc[cortex_ids_int][thal_ids_str] # this gives cx to thal connecitions shape (43, 40)\n",
    "W_thcx_ = full_matrix.loc[thal_ids_int][cortex_ids_str] # this gives thal to cx connecitions shape (40, 43)\n",
    "W = np.array(W_).transpose() # ake a transpose of the W since what the paper used in the transponsed version.\n",
    "W_cxth = np.array(W_cxth_).transpose()\n",
    "W_thcx = np.array(W_thcx_).transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "express-discipline",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:33.203292Z",
     "start_time": "2022-04-20T06:03:33.198874Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# not include RT\n",
    "print(thal_acr[34])\n",
    "thal_acr.remove('RT')\n",
    "plt.imshow(W_cxth)\n",
    "W_cxth = np.delete(W_cxth, 34, 0)\n",
    "print(W_cxth.shape)\n",
    "W_thcx = np.delete(W_thcx, 34, 1)\n",
    "print(W_thcx.shape)\n",
    "\n",
    "# save thalamic arealist\n",
    "print(len(thal_acr))\n",
    "np.savetxt(\"results/thal_list.csv\", thal_acr, delimiter=\",\", fmt = '%s')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "infinite-token",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:34.806137Z",
     "start_time": "2022-04-20T06:03:33.597709Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# peek on the connectivity\n",
    "plot_connectivity(W, cortex_acr,savefig=True)\n",
    "plt.figure()\n",
    "plt.imshow(W>1e-4)\n",
    "W>0 / (len(W)**2 - len(W))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "considered-saskatchewan",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:34.815292Z",
     "start_time": "2022-04-20T06:03:34.808034Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "    # areaidxlist = np.arange(len(cortex_acr))\n",
    "    # # extract the connectivity for list of areas.\n",
    "    # W_ = W[areaidxlist, :][:, areaidxlist]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "victorian-concept",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:34.955743Z",
     "start_time": "2022-04-20T06:03:34.817657Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.imshow(W)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "prompt-edition",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:34.962016Z",
     "start_time": "2022-04-20T06:03:34.958039Z"
    },
    "code_folding": [],
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# generate connectivity matrix \n",
    "# this function also sorted the cortical areas in W_cxth and W_thcx.\n",
    "conn_cxcx,conn_cxth, conn_thcx, area_list_sort, thal_list_sort = generate_connectivity(\n",
    "    W, cortex_acr, W_cxth, W_thcx, thal_acr ,imgout=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "vocal-establishment",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:35.426469Z",
     "start_time": "2022-04-20T06:03:35.423937Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(area_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "charitable-princess",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:35.838955Z",
     "start_time": "2022-04-20T06:03:35.836675Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# # quick check which thalamic area project to cortical area  # SSp-n 26  MOs 14 \n",
    "# np.set_printoptions(precision=3, suppress=True)\n",
    "# print(conn_thcx[14, :]*1e5) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "patent-minutes",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:36.176200Z",
     "start_time": "2022-04-20T06:03:36.158845Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# load the volume data and make processing.\n",
    "dataFile = 'CCFv3SuppTable2.csv'\n",
    "df = pd.read_csv(dataFile, header=1)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "appreciated-wireless",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-01-12T04:33:13.440808Z",
     "start_time": "2022-01-12T04:33:13.436889Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "positive-tutorial",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:36.824642Z",
     "start_time": "2022-04-20T06:03:36.774379Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "volumeList = []\n",
    "for area in area_list:\n",
    "    volume = float(df[df['abbreviation'] == area]['total_voxel_counts (10 um)'])\n",
    "    volumeList.append(volume)\n",
    "volume_df = pd.DataFrame([x for x in zip(area_list, volumeList)], columns=['area', 'volume'])\n",
    "print(volume_df['volume'].max()/volume_df['volume'].min())\n",
    "volume_df.head()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f646d62d",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "smoking-still",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:38.486310Z",
     "start_time": "2022-04-20T06:03:37.171078Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = px.line(volume_df, 'area', 'volume')\n",
    "fig.show()\n",
    "fig.write_image(\"figure/loadDataProcessing/areaVolume.svg\", width=1000, height=350, scale=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dried-scientist",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:38.489846Z",
     "start_time": "2022-04-20T06:03:38.488192Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# load layer 23 cell density\n",
    "# rawl23density_df = pd.read_csv('results/l23cell_density_iso.csv')\n",
    "# rawl23density_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accompanied-number",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:38.494223Z",
     "start_time": "2022-04-20T06:03:38.491977Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# # generate connectivity that is normalized by volume   and using layer 23 cell density\n",
    "# areaV3toV2 = {'VISa':'PTLp', 'FRP':'PL', 'VISli':'VISl', \n",
    "#         'VISrl':'PTLp', 'VISpor':'VISpl', 'SSp-un': 'SSp-ul'}  # copied from helper.py \n",
    "# l23densityList = []\n",
    "# for area in area_list: \n",
    "# #     if curRow.empty: \n",
    "#     if area in areaV3toV2: \n",
    "#         curRow = rawl23density_df.loc[rawl23density_df['Acronym'] == areaV3toV2[area]] \n",
    "#     else:\n",
    "#         curRow = rawl23density_df.loc[rawl23density_df['Acronym'] == area] \n",
    "#     l23densityList.append(float(curRow['Neurons [mm-3]'])) \n",
    "# l23density_df = pd.DataFrame([x for x in zip(area_list, l23densityList)], \n",
    "#                           columns=['area', 'l23 neuron density'])\n",
    "\n",
    "# l23density_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "worldwide-assistant",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:38.509683Z",
     "start_time": "2022-04-20T06:03:38.495829Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# load cell density  paper: A cell atlas for the mouse brain.\n",
    "rawdensityMarkram_df = pd.read_csv('results/cell_density_iso_Markram.csv')\n",
    "rawdensityMarkram_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "collectible-hammer",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:38.552933Z",
     "start_time": "2022-04-20T06:03:38.511500Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#generate connectivity that is normalized by volume and norm by cell density\n",
    "areaV3toV2 = {'VISa':'PTLp', 'FRP':'PL', 'VISli':'VISl', \n",
    "        'VISrl':'PTLp', 'VISpor':'VISpl', 'SSp-un': 'SSp-ul'}  # copied from helper.py \n",
    "densityMarkramList = []\n",
    "exdensityMarkramList = []\n",
    "inhdensityMarkramList = []\n",
    "for area in area_list: \n",
    "#     if curRow.empty: \n",
    "    if area in areaV3toV2: \n",
    "        curRow = rawdensityMarkram_df.loc[rawdensityMarkram_df['Acronym'] == areaV3toV2[area]] \n",
    "    else:\n",
    "        curRow = rawdensityMarkram_df.loc[rawdensityMarkram_df['Acronym'] == area] \n",
    "    densityMarkramList.append(float(curRow['Neurons'])) \n",
    "    exdensityMarkramList.append(float(curRow['Excitatory']))\n",
    "    inhdensityMarkramList.append(float(curRow['Inhibitory']))\n",
    "densityMarkram_df = pd.DataFrame([x for x in zip(area_list, densityMarkramList, exdensityMarkramList, inhdensityMarkramList)], \n",
    "                          columns=['area', 'neuron density', 'excitatory neuron density', 'inhibitory neuron density']) \n",
    "\n",
    "densityMarkram_df.head() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "major-peoples",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "earlier-community",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:38.860226Z",
     "start_time": "2022-04-20T06:03:38.731500Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "fig = px.line(densityMarkram_df, 'area', 'neuron density')\n",
    "fig.show()\n",
    "fig.write_image(\"figure/loadDataProcessing/neuronDensity.svg\", width=1000, height=350, scale=1)\n",
    "\n",
    "densityMarkram_df['ratio'] = densityMarkram_df['excitatory neuron density'] / densityMarkram_df['inhibitory neuron density']\n",
    "fig = px.line(densityMarkram_df, 'area', ['ratio'])\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "durable-knife",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:39.860462Z",
     "start_time": "2022-04-20T06:03:39.841650Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# seems not used. \n",
    "# load cell density\n",
    "newrawdensity_df = pd.read_csv('results/newrawcell_density_iso.csv')\n",
    "newrawdensity_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "flexible-smile",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:42.089844Z",
     "start_time": "2022-04-20T06:03:42.042450Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# This is also from paper: A cell atlas for the mouse brain.   \n",
    "rawdensity_df = pd.read_csv('results/cell_density_iso.csv')  # cell_density_iso.csv is same as cell_density_iso_Markram.csv\n",
    "rawdensity_df.head()\n",
    "\n",
    "areaV3toV2 = {'VISa':'PTLp', 'FRP':'PL', 'VISli':'VISl', \n",
    "        'VISrl':'PTLp', 'VISpor':'VISpl', 'SSp-un': 'SSp-ul'}  # copied from helper.py \n",
    "densityList = []\n",
    "\n",
    "for area in area_list: \n",
    "#     if curRow.empty: \n",
    "    if area in areaV3toV2: \n",
    "#         print(rawdensity_df[rawdensity_df['Acronym'] == areaV3toV2[area]])\n",
    "        curRow = rawdensity_df.loc[rawdensity_df['Acronym'] == areaV3toV2[area]]   \n",
    "#         curRow = rawdensity_df.loc[]\n",
    "    else:\n",
    "        curRow = rawdensity_df.loc[rawdensity_df['Acronym'] == area] \n",
    "    densityList.append(float(curRow['Neurons'])) \n",
    "density_df = pd.DataFrame([x for x in zip(area_list, densityList)], \n",
    "                          columns=['area', 'neuron density']) \n",
    "\n",
    "density_df.head() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "advised-pixel",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:42.578503Z",
     "start_time": "2022-04-20T06:03:42.520301Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = px.line(density_df, 'area', 'neuron density')\n",
    "fig.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sunset-brooks",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:43.240700Z",
     "start_time": "2022-04-20T06:03:43.238925Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# # plot layer 23 density\n",
    "# fig = px.line(l23density_df, 'area', 'l23 neuron density')\n",
    "# fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "current-stick",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:43.802134Z",
     "start_time": "2022-04-20T06:03:43.797864Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# modify the connectivity\n",
    "# mm11\n",
    "normDensity = densityList / np.mean(densityList)\n",
    "densityMat = np.matlib.repmat(normDensity, len(normDensity), 1)\n",
    "\n",
    "# normVolume = np.array(volume_df['volume'] / volume_df['volume'].mean())\n",
    "normVolume = volumeList / np.mean(volumeList)\n",
    "volumeMat = np.matlib.repmat(normVolume, len(normVolume), 1)\n",
    "\n",
    "# conn_cxcx_densityNorm = conn_cxcx / (densityMat * densityMat.T)\n",
    "# conn_cxcx_densityNorm = conn_cxcx / densityMat.T   # just for normalize based on target area density.\n",
    "# conn_cxcx_densityNorm = conn_cxcx * volumeMat / densityMat.T   # normalize based on target area density and by volume\n",
    "conn_cxcx_norm = conn_cxcx * volumeMat / densityMat.T # scaling based on volume\n",
    "\n",
    "# conn_cxcx = conn_cxcx_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "216165b0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "extreme-commissioner",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-19T03:31:38.227796Z",
     "start_time": "2022-04-19T03:31:38.073840Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "appointed-newark",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:46.509109Z",
     "start_time": "2022-04-20T06:03:46.300521Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# sum(conn_cxcx[0, :] * hierarchy.T.squeeze()) / sum(conn_cxcx[0, :]), hierarchy[1][0]\n",
    "plt.imshow(conn_cxcx_norm)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bottom-istanbul",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:47.169284Z",
     "start_time": "2022-04-20T06:03:47.067696Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# get symmetrical connectivity\n",
    "sym_cxcx = 0.5 * (conn_cxcx + conn_cxcx.T)\n",
    "anti_cxcx = 0.5 * (conn_cxcx - conn_cxcx.T)\n",
    "sym = (norm(sym_cxcx) - norm(anti_cxcx)) / (norm(sym_cxcx) + norm(anti_cxcx))\n",
    "print(sym) # if sym ==1 or -1, then conn_cxcx is symmetric. \n",
    "px.imshow(sym_cxcx, color_continuous_scale='viridis')\n",
    "# conn_cxcx_norm\n",
    "sym_cxcx_norm = 0.5 * (conn_cxcx_norm + conn_cxcx_norm.T)\n",
    "px.imshow(sym_cxcx_norm, color_continuous_scale='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "weighted-joyce",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mathematical-church",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:47.969043Z",
     "start_time": "2022-04-20T06:03:47.962787Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "noHierList = [8,9, 11,18,29,32]\n",
    "print([area_list[i] for i in noHierList])\n",
    "h = harrisHierarchy[:, 0]\n",
    "print(h)\n",
    "hierarchy = []\n",
    "j = 0\n",
    "for i in range(n_areas):\n",
    "    if i in noHierList:\n",
    "        hierarchy.append(0)\n",
    "    else:\n",
    "        hierarchy.append(h[j])\n",
    "        j += 1\n",
    "hierarchy = np.array(hierarchy)\n",
    "print(hierarchy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "upper-choir",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:48.300296Z",
     "start_time": "2022-04-20T06:03:48.205504Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sym_cxcx = 0.5 * (conn_cxcx + conn_cxcx.T)\n",
    "anti_cxcx = 0.5 * (conn_cxcx - conn_cxcx.T)\n",
    "sym = (norm(sym_cxcx) - norm(anti_cxcx)) / (norm(sym_cxcx) + norm(anti_cxcx))\n",
    "print(sym) # if sym ==1 or -1, then conn_cxcx is symmetric. \n",
    "px.imshow(sym_cxcx, color_continuous_scale='viridis')\n",
    "# conn_cxcx_norm\n",
    "sym_cxcx_norm = 0.5 * (conn_cxcx_norm + conn_cxcx_norm.T)\n",
    "px.imshow(sym_cxcx_norm, color_continuous_scale='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "colored-pattern",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:49.212852Z",
     "start_time": "2022-04-20T06:03:48.464224Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "conn_withHier = np.zeros((n_areas, n_areas))\n",
    "conn_withHier_norm = np.zeros((n_areas, n_areas))\n",
    "for i in range(n_areas):\n",
    "    for j in range(n_areas):\n",
    "        if j in noHierList:\n",
    "            continue\n",
    "        conn_withHier[i, j] = sym_cxcx[i, j]\n",
    "#         conn_withHier[i, j] = conn_cxcx[j, i]\n",
    "        conn_withHier_norm[i, j] = sym_cxcx_norm[i, j]\n",
    "#         conn_withHier_norm[i, j] = conn_cxcx_norm[j, i]\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(conn_withHier)\n",
    "plt.colorbar()\n",
    "plt.xticks(list(range(n_areas)),area_list,rotation=90,fontname='Georgia',fontsize=7)\n",
    "plt.yticks(list(range(n_areas)),area_list,fontname='Georgia',fontsize=7)\n",
    "plt.show()\n",
    "# conn_6_thres = conn_6 * (conn_6 > threshold_6)\n",
    "# conn_6_norm_thres = conn_6_norm * (conn_6_norm > threshold_6_norm)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "important-facial",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:49.219931Z",
     "start_time": "2022-04-20T06:03:49.214871Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def getPredHi(areaNo, conn, beta):  # not used anymore\n",
    "    s1, s2 = 0, 0\n",
    "    for i in range(n_areas):\n",
    "        if i in noHierList:\n",
    "            continue\n",
    "        s1 += np.exp(beta * conn[areaNo, i]) * hierarchy[i]\n",
    "        s2 += np.exp(beta * conn[areaNo, i])\n",
    "    s = s1 / s2\n",
    "#     h = 0.67 * s + 0.15\n",
    "    h = 0.73 * s + 0.11\n",
    "    return h\n",
    "\n",
    "def getPredHi2(areaNo, conn): \n",
    "    s1, s2 = 0, 0\n",
    "    for i in range(n_areas):\n",
    "        if i in noHierList:\n",
    "            continue\n",
    "        s1 += conn[areaNo, i] * hierarchy[i]\n",
    "        s2 += conn[areaNo, i]\n",
    "    s = s1 / s2\n",
    "    h = 1.33 * s - 0.22\n",
    "    return h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "special-cambodia",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:49.227736Z",
     "start_time": "2022-04-20T06:03:49.222263Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import copy\n",
    "hierarchyEXP = copy.copy(hierarchy)\n",
    "for i in noHierList:\n",
    "    hierarchyEXP[i] = getPredHi(i, conn_cxcx, beta=0.7e4)\n",
    "    hierarchy[i] = getPredHi2(i, conn_cxcx)\n",
    "    print(hierarchy[i], area_list[i])\n",
    "hierarchy_df = pd.DataFrame(hierarchy, index=area_list, columns=['hierarchy index'])\n",
    "# print(hierarchy_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "verified-purple",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:50.427253Z",
     "start_time": "2022-04-20T06:03:49.871023Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "plotHiernotebook(hierarchy_df, 8, r'Hierarchy', (5,3), 12, True, 'loadDataProcessing/hierarchyPlotnotebook.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "separate-ballot",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:03:55.129544Z",
     "start_time": "2022-04-20T06:03:50.429188Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import brian2\n",
    "import pickle\n",
    "import networkx as nx\n",
    "import copy\n",
    "import seaborn as sns\n",
    "from plotFunc import plotPaperTag, getConnHist, getNnodeLoop, plotXYcompr, plotXYcomprCorr, plotXYcomprResidual\n",
    "\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import classification_report, confusion_matrix\n",
    "\n",
    "from adjustText import adjust_text\n",
    "# import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "impressed-graphic",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:18.843462Z",
     "start_time": "2022-04-20T06:03:55.131621Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5, 3),dpi=300)\n",
    "plt.scatter(hierarchy, hierarchyEXP)\n",
    "plt.xlabel('hierarchy(linear)')\n",
    "plt.ylabel('hierarchy(softmax)')\n",
    "fontSize = 12\n",
    "ax = plt.gca()\n",
    "texts = []\n",
    "for i in range(len(area_list)):\n",
    "    acr = area_list[i]\n",
    "    for div_name, div_color in zip(div_name_list, div_color_list):\n",
    "        if acr in div[div_name]:\n",
    "            texts += [ax.text(hierarchy[i], hierarchyEXP[i], acr,\n",
    "                              color=div_color, fontsize=fontSize*0.4)]\n",
    "\n",
    "# # use adjust library to adjust the position of annotations.\n",
    "if True:\n",
    "    adjust_text(texts, hierarchy, hierarchyEXP,\n",
    "                ax=ax, precision=0.001,\n",
    "                arrowprops=dict(arrowstyle='-', color='gray', alpha=.8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "average-california",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-19T04:17:40.941168Z",
     "start_time": "2022-04-19T04:17:40.937514Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "downtown-genius",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:18.848045Z",
     "start_time": "2022-04-20T06:04:18.845316Z"
    },
    "code_folding": [
     0
    ],
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# from scipy.stats import linregress\n",
    "# # linregress(a, b)\n",
    "# i = 8\n",
    "# print(area_list[i])\n",
    "# plt.plot(conn_cxcx[i, :]) # input\n",
    "# plt.figure()\n",
    "# plt.plot(conn_cxcx[:, i]) # output\n",
    "# np.exp(1e4*conn_cxcx[i, :]) / sum(np.exp(1e4*conn_cxcx[i, :]))\n",
    "# np.power(conn_cxcx[i, :], 1/2)\n",
    "\n",
    "# allhierPred, allhierPred2 = [], []\n",
    "# for i in range(n_areas):\n",
    "# #     hPred = sum(np.power(conn_cxcx[i, :], 3) * hierarchy.T.squeeze()) / sum(np.power(conn_cxcx[i, :], 3))\n",
    "#     hPred = sum(np.exp(1e4 * conn_cxcx[i, :]) * hierarchy.T.squeeze()) / sum(np.exp(1e4 * conn_cxcx[i, :]))\n",
    "# #     hPred = sum(conn_cxcx[i, :] * hierarchy.T.squeeze()) / sum(conn_cxcx[i, :])\n",
    "#     hPred2 = sum(conn_cxcx[:, i]**1 * hierarchy.T.squeeze()) / sum(conn_cxcx[:, i]**1)\n",
    "#     allhierPred.append(hPred)\n",
    "#     allhierPred2.append(hPred2)\n",
    "    \n",
    "# allhierPred = np.array(allhierPred)\n",
    "# allhierPred2 = np.array(allhierPred2)\n",
    "# Yhier = hierarchy.T.squeeze()\n",
    "\n",
    "# plt.figure()\n",
    "# plt.scatter(allhierPred, Yhier)\n",
    "# r1, p1 = sp.stats.pearsonr(allhierPred, Yhier)\n",
    "# lg = linregress(allhierPred, Yhier)\n",
    "# plt.title(str(r1))\n",
    "# plt.xlabel('prediction hi')\n",
    "# plt.ylabel('real hi')\n",
    "# plt.xlim(0, 1)\n",
    "\n",
    "\n",
    "# plt.figure()\n",
    "# plt.scatter(allhierPred2, Yhier)\n",
    "# r, p = sp.stats.pearsonr(allhierPred2, Yhier)\n",
    "# plt.title(str(r))\n",
    "# plt.xlabel('prediction hi')\n",
    "# plt.ylabel('real hi')\n",
    "# plt.xlim(0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accompanied-requirement",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:18.851529Z",
     "start_time": "2022-04-20T06:04:18.849690Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# make plot using densityNorm connectivity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mighty-underground",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:18.953972Z",
     "start_time": "2022-04-20T06:04:18.853532Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "sym_cxcx = 0.5 * (conn_cxcx + conn_cxcx.T)\n",
    "anti_cxcx = 0.5 * (conn_cxcx - conn_cxcx.T)\n",
    "sym = (norm(sym_cxcx) - norm(anti_cxcx)) / (norm(sym_cxcx) + norm(anti_cxcx))\n",
    "print(sym) # if sym ==1 or -1, then conn_cxcx is symmetric. \n",
    "px.imshow(sym_cxcx, color_continuous_scale='viridis')\n",
    "# conn_cxcx_norm\n",
    "sym_cxcx_norm = 0.5 * (conn_cxcx_norm + conn_cxcx_norm.T)\n",
    "px.imshow(sym_cxcx_norm, color_continuous_scale='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "massive-cleveland",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:20.311886Z",
     "start_time": "2022-04-20T06:04:18.957612Z"
    },
    "collapsed": true,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(sym_cxcx[sym_cxcx != 0], bins=np.logspace(np.log10(1e-7),np.log10(1e-2), 50), log=True)\n",
    "ax = plt.gca()\n",
    "ax.set_xscale('log')\n",
    "# plt.figure()\n",
    "# plt.hist(conn_cxcx[conn_cxcx != 0], bins=np.logspace(np.log10(1e-7),np.log10(1e-2), 50), log=True)\n",
    "# ax = plt.gca()\n",
    "# ax.set_xscale('log')\n",
    "# plt.figure()\n",
    "# plt.hist(conn_cxcx_norm[conn_cxcx_norm != 0], \n",
    "#          bins=np.logspace(np.log10(1e-7),np.log10(1e-2), 50), log=True)\n",
    "# ax = plt.gca()\n",
    "# ax.set_xscale('log')\n",
    "plt.figure()\n",
    "plt.hist(sym_cxcx_norm[sym_cxcx_norm != 0], \n",
    "         bins=np.logspace(np.log10(1e-7),np.log10(1e-2), 50), log=True)\n",
    "ax = plt.gca()\n",
    "ax.set_xscale('log')\n",
    "\n",
    "\n",
    "# use the percentile\n",
    "\n",
    "print(np.percentile(sym_cxcx, 70))\n",
    "print(np.percentile(sym_cxcx_norm, 65))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "piano-slovak",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:20.319730Z",
     "start_time": "2022-04-20T06:04:20.314496Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# generate a graph for the connectivity\n",
    "# https://github.com/rweng18/midsummer_network\n",
    "conn_graph = nx.Graph()\n",
    "for area in area_list_sort:\n",
    "    conn_graph.add_node(area, size = 2)\n",
    "conn_graph.nodes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "active-sacrifice",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:20.327519Z",
     "start_time": "2022-04-20T06:04:20.322683Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# thresholding the network\n",
    "threshold = 3.21e-4\n",
    "conn_cxcx_thres = sym_cxcx * (sym_cxcx > threshold)\n",
    "# add edges to the graph\n",
    "for area1_index in range(len(area_list_sort)):\n",
    "    for area2_index in range(len(area_list_sort)):\n",
    "        if area1_index > area2_index:\n",
    "            continue\n",
    "        area1 = area_list_sort[area1_index]\n",
    "        area2 = area_list_sort[area2_index]\n",
    "        conn = conn_cxcx_thres[area1_index, area2_index]\n",
    "        if conn > 0:\n",
    "            conn_graph.add_edge(\n",
    "            area1, area2, weight = conn)\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "wireless-contrary",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-11-01T05:32:15.599153Z",
     "start_time": "2021-11-01T05:32:15.594771Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "exempt-refund",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:20.378309Z",
     "start_time": "2022-04-20T06:04:20.329379Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# generate position for each node\n",
    "\n",
    "pos_ = nx.spring_layout(conn_graph)\n",
    "# \n",
    "brainCord = pd.read_csv('brainmapCord.csv')\n",
    "\n",
    "c1 = np.cos(-0.95)\n",
    "s1 = np.sin(-0.95)\n",
    "\n",
    "for area in pos_:\n",
    "    targetRow = brainCord[brainCord['area'] == area]\n",
    "# print(targetRow)\n",
    "    xcord = np.float(targetRow['xcord'])\n",
    "    ycord = np.float(targetRow['ycord'])\n",
    "    xcord_ = xcord * c1 - ycord * s1\n",
    "    ycord_ = xcord * s1 + ycord * c1\n",
    "    pos_[area] = np.array([xcord_, ycord_])\n",
    "\n",
    "# pos_['AIp'] = [16, 5]\n",
    "pos_['MOp'][1] = -5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "competent-format",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:20.383874Z",
     "start_time": "2022-04-20T06:04:20.380101Z"
    },
    "code_folding": [
     0
    ],
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def make_edge(x, y, text, width):\n",
    "    \n",
    "    '''Creates a scatter trace for the edge between x's and y's with given width\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x    : a tuple of the endpoints' x-coordinates in the form, tuple([x0, x1, None])\n",
    "    \n",
    "    y    : a tuple of the endpoints' y-coordinates in the form, tuple([y0, y1, None])\n",
    "    \n",
    "    width: the width of the line\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    An edge trace that goes between x0 and x1 with specified width.\n",
    "    '''\n",
    "    return  go.Scatter(x         = x,\n",
    "                       y         = y,\n",
    "                       line      = dict(width = width,\n",
    "                                   color = 'cornflowerblue'),\n",
    "                       opacity = 0.5,\n",
    "                       hoverinfo = 'text',\n",
    "                       text      = ([text]),\n",
    "                       mode      = 'lines')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "assumed-superintendent",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:20.388466Z",
     "start_time": "2022-04-20T06:04:20.386018Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def weight2thick(wgt, threshold):\n",
    "    return np.log10(wgt) - np.log10(threshold)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "serial-platinum",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:20.458177Z",
     "start_time": "2022-04-20T06:04:20.390482Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# For each edge, make an edge_trace, append to list\n",
    "\n",
    "edge_trace = []\n",
    "for edge in conn_graph.edges():\n",
    "    wgt = conn_graph.edges()[edge]['weight']\n",
    "    if wgt > 0:\n",
    "        area1 = edge[0]\n",
    "        area2 = edge[1]\n",
    "\n",
    "        x0, y0 = pos_[area1]\n",
    "        x1, y1 = pos_[area2]\n",
    "\n",
    "        text   = area1 + '--' + area2 + ': ' + str(wgt)\n",
    "        \n",
    "        trace  = make_edge(\n",
    "            [x0, x1, None], [y0, y1, None], text,\n",
    "            weight2thick(wgt, threshold) * 2 + 0.5)\n",
    "        edge_trace.append(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "focal-scope",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:20.522710Z",
     "start_time": "2022-04-20T06:04:20.459883Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    " # Make a node trace\n",
    "node_trace = go.Scatter(x         = [],\n",
    "                        y         = [],\n",
    "                        text      = [],\n",
    "                        textposition = \"top center\",\n",
    "                        textfont_size = 5,\n",
    "                        mode      = 'markers+text',\n",
    "                        hoverinfo = 'none',\n",
    "                        marker    = dict(color = [],\n",
    "                                         size  = [],\n",
    "                                         line  = None))\n",
    "# For each node, get the position and size and add to the node_trace\n",
    "for node in conn_graph.nodes():\n",
    "    x, y = pos_[node]\n",
    "    node_trace['x'] += tuple([x])\n",
    "    node_trace['y'] += tuple([y])\n",
    "    node_trace['marker']['color'] += tuple(['cornflowerblue'])\n",
    "    node_trace['marker']['size'] += tuple([5*conn_graph.nodes()[node]['size']])\n",
    "    node_trace['text'] += tuple(['<b>' + node + '</b>'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dbf0cde",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the range for thresholding the network.\n",
    "a  = sym_cxcx.reshape(-1,1)\n",
    "plt.hist(a, bins = 100)\n",
    "print(sum(a > 4.2e-4)/len(a)) \n",
    "\n",
    "\n",
    "conn_graph = nx.Graph()\n",
    "for area in area_list_sort:\n",
    "    conn_graph.add_node(area, size = 2)\n",
    "conn_graph.nodes()\n",
    "\n",
    "# thresholding the network\n",
    "threshold = 2.9e-4\n",
    "conn_cxcx_thres = sym_cxcx * (sym_cxcx > threshold)\n",
    "# add edges to the graph\n",
    "for area1_index in range(len(area_list_sort)):\n",
    "    for area2_index in range(area1_index+1, len(area_list_sort)):\n",
    "        area1 = area_list_sort[area1_index]\n",
    "        area2 = area_list_sort[area2_index]\n",
    "        conn = conn_cxcx_thres[area1_index, area2_index]\n",
    "        if conn > 0:\n",
    "            conn_graph.add_edge(\n",
    "            area1, area2, weight = conn)\n",
    "\n",
    "# For each edge, make an edge_trace, append to list\n",
    "\n",
    "edge_trace = []\n",
    "for edge in conn_graph.edges():\n",
    "    wgt = conn_graph.edges()[edge]['weight']\n",
    "    if wgt > 0:\n",
    "        area1 = edge[0]\n",
    "        area2 = edge[1]\n",
    "\n",
    "        x0, y0 = pos_[area1]\n",
    "        x1, y1 = pos_[area2]\n",
    "\n",
    "        text   = area1 + '--' + area2 + ': ' + str(wgt)\n",
    "        \n",
    "        trace  = make_edge(\n",
    "            [x0, x1, None], [y0, y1, None], text,\n",
    "            weight2thick(wgt, threshold) * 2 + 0.5)\n",
    "        edge_trace.append(trace)\n",
    "        \n",
    "# Make a node trace\n",
    "node_trace = go.Scatter(x         = [],\n",
    "                        y         = [],\n",
    "                        text      = [],\n",
    "                        textposition = \"top center\",\n",
    "                        textfont_size = 5,\n",
    "                        mode      = 'markers+text',\n",
    "                        hoverinfo = 'none',\n",
    "                        marker    = dict(color = [],\n",
    "                                         size  = [],\n",
    "                                         line  = None))\n",
    "# For each node, get the position and size and add to the node_trace\n",
    "for node in conn_graph.nodes():\n",
    "    x, y = pos_[node]\n",
    "    node_trace['x'] += tuple([x])\n",
    "    node_trace['y'] += tuple([y])\n",
    "    node_trace['marker']['color'] += tuple(['cornflowerblue'])\n",
    "    node_trace['marker']['size'] += tuple([5*conn_graph.nodes()[node]['size']])\n",
    "    node_trace['text'] += tuple(['<b>' + node + '</b>'])        \n",
    "\n",
    "\n",
    "layout = go.Layout(\n",
    "    paper_bgcolor='rgba(0,0,0,0)',\n",
    "    plot_bgcolor='rgba(0,0,0,0)'\n",
    ")\n",
    "\n",
    "saveFig = True\n",
    "fileName = 'rawConnectivityMap.svg'\n",
    "\n",
    "fig = go.Figure(layout = layout)\n",
    "\n",
    "for trace in edge_trace:\n",
    "    fig.add_trace(trace)\n",
    "\n",
    "fig.add_trace(node_trace)\n",
    "\n",
    "fig.update_layout(showlegend = False)\n",
    "\n",
    "fig.update_xaxes(showticklabels = False, automargin=True)\n",
    "\n",
    "fig.update_yaxes(showticklabels = False, scaleanchor='x', scaleratio=1, automargin=True)\n",
    "fig.update_layout(width = 700, height=580)\n",
    "fig.show()\n",
    "if saveFig:\n",
    "    fig.write_image(\"figure/loadDataProcessing/\"+ fileName, width=700, height=580, scale=1)\n",
    "# py.plot(fig, filename='conn_graph.html')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "romantic-initial",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:21.771835Z",
     "start_time": "2022-04-20T06:04:20.524346Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "conn_graph = nx.Graph()\n",
    "for area in area_list_sort:\n",
    "    conn_graph.add_node(area, size = 2)\n",
    "conn_graph.nodes()\n",
    "\n",
    "# thresholding the network\n",
    "threshold = 4.2e-4 #2.9e-4\n",
    "conn_cxcx_thres = sym_cxcx * (sym_cxcx > threshold)\n",
    "# add edges to the graph\n",
    "for area1_index in range(len(area_list_sort)):\n",
    "    for area2_index in range(area1_index+1, len(area_list_sort)):\n",
    "        area1 = area_list_sort[area1_index]\n",
    "        area2 = area_list_sort[area2_index]\n",
    "        conn = conn_cxcx_thres[area1_index, area2_index]\n",
    "        if conn > 0:\n",
    "            conn_graph.add_edge(\n",
    "            area1, area2, weight = conn)\n",
    "\n",
    "# For each edge, make an edge_trace, append to list\n",
    "\n",
    "edge_trace = []\n",
    "for edge in conn_graph.edges():\n",
    "    wgt = conn_graph.edges()[edge]['weight']\n",
    "    if wgt > 0:\n",
    "        area1 = edge[0]\n",
    "        area2 = edge[1]\n",
    "\n",
    "        x0, y0 = pos_[area1]\n",
    "        x1, y1 = pos_[area2]\n",
    "\n",
    "        text   = area1 + '--' + area2 + ': ' + str(wgt)\n",
    "        \n",
    "        trace  = make_edge(\n",
    "            [x0, x1, None], [y0, y1, None], text,\n",
    "            weight2thick(wgt, threshold) * 2 + 0.5)\n",
    "        edge_trace.append(trace)\n",
    "        \n",
    "# Make a node trace\n",
    "node_trace = go.Scatter(x         = [],\n",
    "                        y         = [],\n",
    "                        text      = [],\n",
    "                        textposition = \"top center\",\n",
    "                        textfont_size = 8,\n",
    "                        mode      = 'markers+text',\n",
    "                        hoverinfo = 'none',\n",
    "                        marker    = dict(color = [],\n",
    "                                         size  = [],\n",
    "                                         line  = None))\n",
    "# For each node, get the position and size and add to the node_trace\n",
    "for node in conn_graph.nodes():\n",
    "    x, y = pos_[node]\n",
    "    node_trace['x'] += tuple([x])\n",
    "    node_trace['y'] += tuple([y])\n",
    "    node_trace['marker']['color'] += tuple(['cornflowerblue'])\n",
    "    node_trace['marker']['size'] += tuple([5*conn_graph.nodes()[node]['size']])\n",
    "    node_trace['text'] += tuple(['<b>' + node + '</b>'])        \n",
    "\n",
    "\n",
    "layout = go.Layout(\n",
    "    paper_bgcolor='rgba(0,0,0,0)',\n",
    "    plot_bgcolor='rgba(0,0,0,0)'\n",
    ")\n",
    "\n",
    "saveFig = True\n",
    "fileName = 'rawConnectivityMap.svg'\n",
    "\n",
    "fig = go.Figure(layout = layout)\n",
    "\n",
    "for trace in edge_trace:\n",
    "    fig.add_trace(trace)\n",
    "\n",
    "fig.add_trace(node_trace)\n",
    "\n",
    "fig.update_layout(showlegend = False)\n",
    "\n",
    "fig.update_xaxes(showticklabels = False, automargin=True)\n",
    "\n",
    "fig.update_yaxes(showticklabels = False, scaleanchor='x', scaleratio=1, automargin=True)\n",
    "fig.update_layout(width = 700, height=580)\n",
    "fig.show()\n",
    "if saveFig:\n",
    "    fig.write_image(\"figure/loadDataProcessing/\"+ fileName, width=700, height=580, scale=1)\n",
    "# py.plot(fig, filename='conn_graph.html')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe486c2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the range for thresholding the network.\n",
    "a  = sym_cxcx_norm.reshape(-1,1)\n",
    "plt.hist(a, bins = 100)\n",
    "print(sum(a > 3.5e-4)/len(a)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "divine-warrant",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:22.420366Z",
     "start_time": "2022-04-20T06:04:21.774130Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# generate a graph for the connectivity\n",
    "# https://github.com/rweng18/midsummer_network\n",
    "\n",
    "saveFig = True\n",
    "fileName = 'connectivityNormalized.svg'\n",
    "\n",
    "conn_graph = nx.Graph()\n",
    "for area in area_list_sort:\n",
    "    conn_graph.add_node(area, size = 2)\n",
    "conn_graph.nodes()\n",
    "\n",
    "# thresholding the network\n",
    "threshold = 3.5e-4  #1.6e-4\n",
    "conn_cxcx_norm_thres = sym_cxcx_norm * (sym_cxcx_norm > threshold)\n",
    "# add edges to the graph\n",
    "for area1_index in range(len(area_list_sort)):\n",
    "    for area2_index in range(area1_index+1, len(area_list_sort)):\n",
    "\n",
    "        area1 = area_list_sort[area1_index]\n",
    "        area2 = area_list_sort[area2_index]\n",
    "        conn = conn_cxcx_norm_thres[area1_index, area2_index]\n",
    "        if conn > 0:\n",
    "            conn_graph.add_edge(\n",
    "            area1, area2, weight = conn)\n",
    "            \n",
    "# For each edge, make an edge_trace, append to list\n",
    "\n",
    "edge_trace = []\n",
    "for edge in conn_graph.edges():\n",
    "    wgt = conn_graph.edges()[edge]['weight']\n",
    "    if wgt > 0:\n",
    "        area1 = edge[0]\n",
    "        area2 = edge[1]\n",
    "\n",
    "        x0, y0 = pos_[area1]\n",
    "        x1, y1 = pos_[area2]\n",
    "\n",
    "        text   = area1 + '--' + area2 + ': ' + str(wgt)\n",
    "        \n",
    "        trace  = make_edge(\n",
    "            [x0, x1, None], [y0, y1, None], text,\n",
    "            weight2thick(wgt, threshold) * 2 + 0.5)\n",
    "        edge_trace.append(trace)\n",
    "        \n",
    "        \n",
    "# Make a node trace\n",
    "node_trace = go.Scatter(x         = [],\n",
    "                        y         = [],\n",
    "                        text      = [],\n",
    "                        textposition = \"top center\",\n",
    "                        textfont_size = 8,\n",
    "                        mode      = 'markers+text',\n",
    "                        hoverinfo = 'none',\n",
    "                        marker    = dict(color = [],\n",
    "                                         size  = [],\n",
    "                                         line  = None))\n",
    "# For each node, get the position and size and add to the node_trace\n",
    "for node in conn_graph.nodes():\n",
    "    x, y = pos_[node]\n",
    "    node_trace['x'] += tuple([x])\n",
    "    node_trace['y'] += tuple([y])\n",
    "    node_trace['marker']['color'] += tuple(['cornflowerblue'])\n",
    "    node_trace['marker']['size'] += tuple([5*conn_graph.nodes()[node]['size']])\n",
    "    node_trace['text'] += tuple(['<b>' + node + '</b>'])   \n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    layout = go.Layout(\n",
    "    paper_bgcolor='rgba(0,0,0,0)',\n",
    "    plot_bgcolor='rgba(0,0,0,0)'\n",
    ")\n",
    "\n",
    "\n",
    "\n",
    "fig = go.Figure(layout = layout)\n",
    "\n",
    "for trace in edge_trace:\n",
    "    fig.add_trace(trace)\n",
    "\n",
    "fig.add_trace(node_trace)\n",
    "\n",
    "fig.update_layout(showlegend = False)\n",
    "\n",
    "fig.update_xaxes(showticklabels = False)\n",
    "\n",
    "fig.update_yaxes(showticklabels = False, scaleanchor='x', scaleratio=1)\n",
    "fig.update_layout(width = 700, height=580)\n",
    "fig.show()\n",
    "if saveFig:\n",
    "    fig.write_image(\"figure/loadDataProcessing/\"+ fileName, width=700, height=580, scale=1)\n",
    "# py.plot(fig, filename='conn_graph.html')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "noted-stability",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:22.519024Z",
     "start_time": "2022-04-20T06:04:22.422213Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "inDegree, outDegree = [], []\n",
    "thresIn, thresOut = 2.75e-4, 2.75e-4\n",
    "for i in range(len(area_list)):\n",
    "    inDegree.append(sum(conn_cxcx_norm[i, :] > thresIn))\n",
    "    outDegree.append(sum(conn_cxcx_norm[:, i] > thresOut))\n",
    "df = pd.DataFrame({'area': area_list, 'inDegree': inDegree, 'outDegree': outDegree})\n",
    "fig = px.line(df, 'area', ['inDegree', 'outDegree'])\n",
    "fig.update_yaxes(title='degree')\n",
    "fig.update_layout(legend_title_text='')\n",
    "fig.show()\n",
    "fig.write_image(\"figure/loadDataProcessing/\"+ \"inOutDegree.svg\", width=1000, height=380, scale=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "horizontal-affect",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:22.524428Z",
     "start_time": "2022-04-20T06:04:22.520707Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# use 75% to determine threshold.\n",
    "print(np.percentile(sym_cxcx_norm, 75))\n",
    "print(np.percentile(sym_cxcx, 75))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6f8612d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeHeatPlotRect(mat, xTicks, yTicks, Xlabel, Ylabel, cbarlabel, colormap, figureSize, fontSize, saveFig, fileName):\n",
    "    fig = plt.figure(figsize=figureSize, dpi=300, facecolor='w', edgecolor='k')\n",
    "    plt.rcParams.update({'font.size': fontSize})\n",
    "#     plt.imshow(mat, colormap, norm=TwoSlopeNorm(0))\n",
    "    plt.imshow(mat, colormap)\n",
    "    cbar = plt.colorbar()\n",
    "    cbar.set_label(cbarlabel)\n",
    "    plt.xticks(range(len(xTicks)), xTicks, rotation=90, fontsize=fontSize*0.4)\n",
    "    plt.yticks(range(len(yTicks)), yTicks, fontsize=fontSize*0.4)\n",
    "    plt.xlabel(Xlabel)\n",
    "    plt.ylabel(Ylabel)\n",
    "    \n",
    "    # save figure\n",
    "    if saveFig:\n",
    "        plt.savefig('figure/' + fileName, dpi=300, bbox_inches='tight', transparent=True)\n",
    "    # plt.savefig('figure/neurondensity_grad.png',dpi=80,bbox_inc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "circular-memory",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:22.657118Z",
     "start_time": "2022-04-20T06:04:22.529063Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# rearrange the connnectivity to be hierarhchical. \n",
    "sortedAreaList = [area_list[i] for i in np.argsort(list(hierarchy_df['hierarchy index']))]\n",
    "print(sortedAreaList[:5], len(sortedAreaList))\n",
    "\n",
    "sortedConnPre = np.array( [conn_cxcx_norm[:,j].tolist() for j in np.argsort(list(hierarchy_df['hierarchy index']))] )\n",
    "sortedConn = np.array( [sortedConnPre[:,i].tolist() for i in np.argsort(list(hierarchy_df['hierarchy index']))] )\n",
    "\n",
    "print(sortedConn)\n",
    "\n",
    "makeHeatPlotRect(sortedConn, sortedAreaList, sortedAreaList,\n",
    "                'Source', 'Target', 'Connectivity', 'viridis', (5, 3), 12, True, 'loadDataProcessing/cortexConnectivityMatrix.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "empirical-bangladesh",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:22.781098Z",
     "start_time": "2022-04-20T06:04:22.659242Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sortedconncxth = np.array( [conn_cxth[:,j].tolist() for j in np.argsort(list(hierarchy_df['hierarchy index']))] ).T\n",
    "# sortedConn = np.array( [sortedConnPre[:,i].tolist() for i in np.argsort(list(hierarchy_df['hierarchy index']))] )\n",
    "# TODO: check the order of the connectivity\n",
    "#  unfinished sorted CXTH connections. \n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(conn_cxth)\n",
    "plt.xticks(range(len(area_list)), area_list, rotation=90)\n",
    "plt.yticks(range(len(thal_acr)), thal_acr)\n",
    "plt.xlabel('source')\n",
    "plt.ylabel('target')\n",
    "plt.colorbar()\n",
    "\n",
    "makeHeatPlotRect(sortedconncxth, sortedAreaList, thal_acr,\n",
    "                'Source', 'Target', 'Connectivity', 'viridis', (5, 3), 12, True, 'loadDataProcessing/cxthConnectivityMatrix.svg')\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5046469",
   "metadata": {},
   "outputs": [],
   "source": [
    "sortedconnthcx = np.array( [conn_thcx[i,:].tolist() for i in np.argsort(list(hierarchy_df['hierarchy index']))] )\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(conn_thcx)\n",
    "plt.yticks(range(len(area_list)), area_list)\n",
    "plt.xticks(range(len(thal_acr)), thal_acr, rotation=90)\n",
    "plt.xlabel('source')\n",
    "plt.ylabel('target')\n",
    "plt.colorbar()\n",
    "\n",
    "makeHeatPlotRect(sortedconnthcx, thal_acr, sortedAreaList,\n",
    "                'Source', 'Target', 'Connectivity', 'viridis', (5, 3), 12, True, 'loadDataProcessing/thcxConnectivityMatrix.svg')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sixth-silly",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:22.788319Z",
     "start_time": "2022-04-20T06:04:22.782892Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# save arealist\n",
    "np.savetxt(\"results/area_list.csv\", sortedAreaList, delimiter=\",\", fmt = '%s')\n",
    "# check if loading area list works\n",
    "arr = np.loadtxt('results/area_list.csv', delimiter=',', dtype='str')\n",
    "arr[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faa16c96",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(sortedConn)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40f6abb8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "packed-hearing",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:22.804206Z",
     "start_time": "2022-04-20T06:04:22.790002Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# save connectivity\n",
    "np.savetxt(\"results/conn_cxcx.csv\", sortedConn, delimiter=\",\")\n",
    "np.savetxt(\"results/conn_cxth.csv\", sortedconncxth, delimiter=\",\")\n",
    "np.savetxt(\"results/conn_thcx.csv\", sortedconnthcx, delimiter=\",\")\n",
    "\n",
    "np.savetxt(\"results/conn_cxcx_raw.csv\", conn_cxcx, delimiter=\",\")\n",
    "\n",
    "\n",
    "volume_df.to_csv('results/dfVolume.csv')\n",
    "density_df.to_csv(\"results/dfDensity.csv\")\n",
    "\n",
    "# np.savetxt(\"results/dfDensity.csv\", density_df, delimiter=\",\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "vulnerable-glossary",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:22.811766Z",
     "start_time": "2022-04-20T06:04:22.805842Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "arr = np.loadtxt('results/conn_cxcx.csv', delimiter=',')\n",
    "print(arr[0:5, 0:5])\n",
    "print(arr.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "macro-evening",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "strategic-booking",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:49.207466Z",
     "start_time": "2022-04-20T06:04:49.194835Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# load interneuron density \n",
    "# and add processing afterwards. \n",
    "\n",
    "interneuronfilename = 'mouse_interneuron_density.pkl'\n",
    "celldensityfilename = 'cell_density_iso.csv'\n",
    "# get density for related old areas (the area name is based on CCF v2, not latest CCF v3)\n",
    "old_areas = ['ACAd', 'ACAv', 'AId', 'AIp', 'AIv', 'AUDd', 'AUDp','AUDpo', 'AUDv', 'ECT', 'GU',\n",
    "'ILA', 'MOp', 'MOs', 'ORBl', 'ORBm', 'ORBvl', 'PERI', 'PL', 'PTLp', 'RSPagl', 'RSPd',\n",
    "'RSPv', 'SSp-bfd', 'SSp-ll', 'SSp-m', 'SSp-n', 'SSp-tr', 'SSp-ul', 'SSs', 'TEa', 'VISC', \n",
    "'VISal', 'VISam', 'VISl', 'VISp', 'VISpl', 'VISpm']\n",
    "\n",
    "neurons_sorted = load_celldensity(libfilepath, celldensityfilename, old_areas)\n",
    " \n",
    "[PVgrad_df, SSTgrad_df, normPVgrad_df, normSSTgrad_df] = load_interneurondensity(\n",
    "    libfilepath, interneuronfilename, old_areas)\n",
    "\n",
    "# modify cell density\n",
    "# neurons_sorted[37] = neurons_sorted[33] # VISpm seems outlier, replace it with VISam\n",
    "\n",
    "# load pandas dataframe \"results/dfDensity.csv\"\n",
    "density_df = pd.read_csv(\"results/dfDensity.csv\", index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ceb5848",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78c8cd49",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eleven-thumb",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:50.364059Z",
     "start_time": "2022-04-20T06:04:50.233746Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "a = np.array(list(PVgrad_df['raw PV density']))\n",
    "r1 = (a - min(a)) / (max(a) - min(a))\n",
    "r2 = a / max(a)\n",
    "plt.scatter(r1, r2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "satisfactory-selection",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:55.392441Z",
     "start_time": "2022-04-20T06:04:55.389747Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "len(neurons_sorted)\n",
    "print(old_areas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "specialized-problem",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:56.441645Z",
     "start_time": "2022-04-20T06:04:56.428202Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "PVcellNorm = np.array(PVgrad_df)[:, 0] / neurons_sorted\n",
    "PVcell_df = pd.DataFrame(PVcellNorm, columns=['raw PV density'], index=old_areas)\n",
    "PVnormcellNorm = np.array(normPVgrad_df)[:, 0] / neurons_sorted \n",
    "PVnormcellNorm = PVnormcellNorm / PVnormcellNorm.max()\n",
    "normPVcell_df = pd.DataFrame(PVnormcellNorm, columns=['norm PV gradient'], index=old_areas)\n",
    "\n",
    "# SSTcell_df is the SST cell fraction. \n",
    "SSTcellNorm = np.array(SSTgrad_df)[:, 0] / neurons_sorted\n",
    "SSTcell_df = pd.DataFrame(SSTcellNorm, columns=['raw SST density'], index=old_areas)\n",
    "SSTcell_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a26806bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "areasInterest = ['SSp-bfd','VISp', 'MOs', 'PL', 'ILA']\n",
    "print([PVcell_df['raw PV density'][A] for A in areasInterest])\n",
    "print([SSTcell_df['raw SST density'][A] for A in areasInterest])\n",
    "\n",
    "# SSp-bfd VIP: 0.071/2 = 0.0355 total: 0.071+0.056+0.0355 = 0.1625\n",
    "# VISp VIP: 0.05*0.6 = 0.03 total: 0.051+0.05+0.03 = 0.131\n",
    "# MOs VIP: 0.056/2 = 0.028   total: 0.056+0.061+0.028 = 0.145\n",
    "# PL VIP: 0.02*1.3 = 0.026 total: 0.02+0.065+0.026 = 0.111\n",
    "# ILA VIP: 0.021*0.8 = 0.0168 total: 0.021+0.089+0.0168 = 0.1268\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "wooden-attraction",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:04:57.918455Z",
     "start_time": "2022-04-20T06:04:57.670929Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(neurons_sorted)\n",
    "plt.xticks(range(len(neurons_sorted)), old_areas, rotation = 90)\n",
    "# plt.plot(sorted(neurons_sorted))\n",
    "# sortIndex = np.argsort(neurons_sorted)\n",
    "# plt.xticks(range(len(neurons_sorted)), [old_areas[i] for i in sortIndex], rotation = 90)\n",
    "plt.ylabel('neuron density')\n",
    "plt.ylim((0, 160000))\n",
    "axes = plt.gca()\n",
    "axes.spines['top'].set_visible(False)\n",
    "axes.spines['right'].set_visible(False)\n",
    "print(sorted(neurons_sorted)[-1] / sorted(neurons_sorted)[0], sorted(neurons_sorted)[-2] / sorted(neurons_sorted)[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "capable-season",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:10.754940Z",
     "start_time": "2022-04-20T06:05:10.692904Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# preprocessing the interneuron density\n",
    "[PVgrad_df, SSTgrad_df, PVcell_df, normPVgrad_df, normSSTgrad_df, normPVcell_df] = interneuronDensityProcessing(\n",
    "    [PVgrad_df, SSTgrad_df, PVcell_df, normPVgrad_df, normSSTgrad_df, normPVcell_df])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dramatic-venice",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:11.251618Z",
     "start_time": "2022-04-20T06:05:11.247607Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(PVcell_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0c8d085",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "velvet-punishment",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:12.277788Z",
     "start_time": "2022-04-20T06:05:12.269097Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def plotRawDensitynotebook(grad_df, dotSize, yaxis, figureSize, fontSize, saveFig, fileName):\n",
    "        # seaborn.set()\n",
    "    dfRaw = grad_df.sort_values(grad_df.columns[0])\n",
    "\n",
    "    fig = plt.figure(figsize = figureSize, dpi=300, facecolor='w', edgecolor='k')\n",
    "    plt.rcParams.update({'font.size': fontSize})\n",
    "    \n",
    "    xx = np.arange(len(grad_df))\n",
    "    plt.plot(xx, dfRaw.values, '.-', markersize=dotSize, color='k')\n",
    "    plt.xticks(xx, dfRaw.index, rotation = 90, fontsize=fontSize*0.6)\n",
    "    plt.ylabel(yaxis) #r'$e_i/e_1$'\n",
    "    # regPV = sk.linear_model.LinearRegression()\n",
    "    # xx2 = xx[np.newaxis].T\n",
    "    # regPV.fit(xx2, dfRaw.values)\n",
    "    # print(regPV.coef_, regPV.intercept_)\n",
    "    # plt.plot(xx, xx * regPV.coef_[0][0] + regPV.intercept_[0], alpha = 0.5)\n",
    "    ax = plt.gca()\n",
    "#     [t.set_color('red') for t in ax.xaxis.get_ticklabels()]\n",
    "    for t in ax.xaxis.get_ticklabels():\n",
    "        label = t._text\n",
    "        for i in range(6):\n",
    "            if label in div[div_name_list[i]]:\n",
    "                t.set_color(div_color_list[i])\n",
    "    # plt.grid(alpha = 0.5)\n",
    "\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    if saveFig:\n",
    "        plt.savefig('figure/' + fileName, dpi=300, bbox_inches='tight',transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac90e1ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotDensityBasedOnRawPVnotebook(grad_df, dotSize, yaxis, figureSize, fontSize, saveFig, fileName):\n",
    "        # seaborn.set()\n",
    "    dfSort = grad_df.sort_values(by = 'raw PV density')\n",
    "    \n",
    "\n",
    "    fig = plt.figure(figsize = figureSize, dpi=300, facecolor='w', edgecolor='k')\n",
    "    plt.rcParams.update({'font.size': fontSize})\n",
    "    \n",
    "    xx = np.arange(len(grad_df))\n",
    "    plt.plot(xx, dfSort['neuron density'], '.-', markersize=dotSize, color='k')\n",
    "    plt.xticks(xx, dfSort['area'], rotation = 90, fontsize=fontSize*0.6)\n",
    "    plt.ylabel(yaxis) #r'$e_i/e_1$'\n",
    "    # regPV = sk.linear_model.LinearRegression()\n",
    "    # xx2 = xx[np.newaxis].T\n",
    "    # regPV.fit(xx2, dfRaw.values)\n",
    "    # print(regPV.coef_, regPV.intercept_)\n",
    "    # plt.plot(xx, xx * regPV.coef_[0][0] + regPV.intercept_[0], alpha = 0.5)\n",
    "    ax = plt.gca()\n",
    "#     [t.set_color('red') for t in ax.xaxis.get_ticklabels()]\n",
    "    for t in ax.xaxis.get_ticklabels():\n",
    "        label = t._text\n",
    "        for i in range(6):\n",
    "            if label in div[div_name_list[i]]:\n",
    "                t.set_color(div_color_list[i])\n",
    "    # plt.grid(alpha = 0.5)\n",
    "\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    if saveFig:\n",
    "        plt.savefig('figure/' + fileName, dpi=300, bbox_inches='tight',transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a12f045",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the cell density according to the order of raw PV density. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alleged-briefs",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:13.412058Z",
     "start_time": "2022-04-20T06:05:12.822852Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "plotRawDensitynotebook(PVgrad_df, 8, r'PV cell density (cells/$mm^{3}$)', (5,3), 12, True, 'loadDataProcessing/rawPVnotebook.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35243591",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acoustic-campaign",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:14.186935Z",
     "start_time": "2022-04-20T06:05:13.638143Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "plotRawDensitynotebook(PVcell_df, 8, r'PV cell fraction', (5,3), 12, True, 'loadDataProcessing/ratioPVnotebook.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "095fa951",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a dataframe for neurondensity in density_df and raw pv density in PVgrad_df.\n",
    "density_rawPV_df = pd.DataFrame({'area': list(density_df['area']), 'neuron density': list(density_df['neuron density']), 'raw PV density': list(PVgrad_df['raw PV density'])})\n",
    "density_rawPV_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2bee7d9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82ff33ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sort the dataframe by raw PV density\n",
    "dfSort = density_rawPV_df.sort_values(by = 'raw PV density')\n",
    "dfSort  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17a9f6a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotDensityBasedOnRawPVnotebook(density_rawPV_df, 8, r'Neuron density (cells/$mm^{3}$)', (5,3), 12, True, 'loadDataProcessing/densityBasedOnPVnotebook.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2f2ad85",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the neuron density with same order as pv cell density.\n",
    "def plotRawDensitynotebook(grad_df, dotSize, yaxis, figureSize, fontSize, saveFig, fileName):\n",
    "    \n",
    "    dfRaw = grad_df.sort_values(grad_df.columns[0])\n",
    "\n",
    "    fig = plt.figure(figsize = figureSize, dpi=300, facecolor='w', edgecolor='k')\n",
    "    plt.rcParams.update({'font.size': fontSize})\n",
    "    \n",
    "    xx = np.arange(len(grad_df))\n",
    "    plt.plot(xx, dfRaw.values, '.-', markersize=dotSize, color='k')\n",
    "    plt.xticks(xx, dfRaw.index, rotation = 90, fontsize=fontSize*0.6)\n",
    "    plt.ylabel(yaxis) #r'$e_i/e_1$'\n",
    "    ax = plt.gca()\n",
    "    for t in ax.xaxis.get_ticklabels():\n",
    "        label = t._text\n",
    "        for i in range(6):\n",
    "            if label in div[div_name_list[i]]:\n",
    "                t.set_color(div_color_list[i])\n",
    "    # plt.grid(alpha = 0.5)\n",
    "\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    if saveFig:\n",
    "        plt.savefig('figure/' + fileName, dpi=300, bbox_inches='tight',transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "homeless-ladder",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:14.779893Z",
     "start_time": "2022-04-20T06:05:14.778025Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# # plot pv norm by l23 density\n",
    "# cellnormPVgrad_df = pd.DataFrame(np.array(PVgrad_df['raw PV density']) / np.array(l23density_df['l23 neuron density']), index = area_list, columns = ['cell norm PV density'])\n",
    "# plotRawDensitynotebook(cellnormPVgrad_df, 8, r'cellnormPV cell density (cells/$mm^{3}$)', (5,3), 12, False, 'loadDataProcessing/cellnormPVnotebook.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hybrid-tomorrow",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:15.251950Z",
     "start_time": "2022-04-20T06:05:15.250032Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# inhcellnormPVgrad_df = pd.DataFrame(np.array(PVgrad_df['raw PV density']) / np.array(density_df['inhibitory neuron density']), index = area_list, columns = ['Inhibitory cell norm PV density'])\n",
    "# plotRawDensitynotebook(inhcellnormPVgrad_df, 8, r'inhibitory cell norm PV density (cells/$mm^{3}$)', (5,3), 12, False, 'loadDataProcessing/inhcellnormPVnotebook.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "unsigned-harmony",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:16.334712Z",
     "start_time": "2022-04-20T06:05:15.765262Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "plotRawDensitynotebook(SSTgrad_df, 8, r'SST cell density (cells/$mm^{3}$)', (5,3), 12, True, 'loadDataProcessing/rawSSTnotebook.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "realistic-appointment",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:16.985574Z",
     "start_time": "2022-04-20T06:05:16.516760Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "dfSST = SSTgrad_df.sort_values(SSTgrad_df.columns[0])\n",
    "\n",
    "plt.figure(figsize = (10, 7))\n",
    "xx = np.arange(len(dfSST))\n",
    "plt.plot(xx, dfSST.values)\n",
    "plt.xticks(xx, dfSST.index, rotation = 90, size = 15)\n",
    "regSST = sk.linear_model.LinearRegression()\n",
    "xx2 = xx[np.newaxis].T\n",
    "regSST.fit(xx2, dfSST.values)\n",
    "print(regSST.coef_, regSST.intercept_)\n",
    "plt.plot(xx, xx * regSST.coef_[0][0] + regSST.intercept_[0], alpha = 0.5)\n",
    "\n",
    "plt.ylabel('SST density (cells/' + r'$mm^{3}$' + ')')\n",
    "plt.grid(alpha = 0.5)\n",
    "axes = plt.gca()\n",
    "axes.spines['top'].set_visible(False)\n",
    "axes.spines['right'].set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "realistic-entity",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:17.271015Z",
     "start_time": "2022-04-20T06:05:17.262252Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "normPVcell_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "located-buffalo",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:44.004076Z",
     "start_time": "2022-04-20T06:05:43.994439Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# generate sorted dataframe of PV density.\n",
    "\n",
    "sortednormPV_df = pd.DataFrame( [normPVcell_df['norm PV gradient'][i] for i in np.argsort(list(hierarchy_df['hierarchy index'])) ], index=sortedAreaList, columns=['norm PV gradient'])\n",
    "sortednormPV_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accompanied-content",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:45.573330Z",
     "start_time": "2022-04-20T06:05:45.565747Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sortedrawPV_df = pd.DataFrame( [PVcell_df['raw PV density'][i] for i in np.argsort(list(hierarchy_df['hierarchy index'])) ], index=sortedAreaList, columns=['raw PV relative density'])\n",
    "sortedrawPV_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "corporate-gender",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:46.686726Z",
     "start_time": "2022-04-20T06:05:46.679507Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# save interneuron density\n",
    "sortednormPV_df.to_csv('results/dfPV.csv')\n",
    "sortedrawPV_df.to_csv('results/dfPVraw.csv')\n",
    "\n",
    "# normPVgrad_df.to_csv('results/dfPV.csv')\n",
    "# normSSTgrad_df.to_csv('results/dfSST.csv')\n",
    "# PVgrad_df.to_csv('results/dfPVraw.csv')\n",
    "# SSTgrad_df.to_csv('results/dfSSTraw.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "understanding-nicholas",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:47.410849Z",
     "start_time": "2022-04-20T06:05:47.400329Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dfPV = pd.read_csv('results/dfPV.csv', header = 0, index_col=0)\n",
    "print(len(dfPV))\n",
    "dfPV.head()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "owned-jungle",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:48.789933Z",
     "start_time": "2022-04-20T06:05:48.780038Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# generate sorted dataframe of hierarchy.\n",
    "\n",
    "sortedHierarchy_df = pd.DataFrame( [hierarchy_df['hierarchy index'][i] for i in np.argsort(list(hierarchy_df['hierarchy index'])) ], index=sortedAreaList, columns=['hierarchy index'])\n",
    "sortedHierarchy_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "center-clothing",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-20T06:05:50.854012Z",
     "start_time": "2022-04-20T06:05:50.849082Z"
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sortedHierarchy_df.to_csv('results/dfHier.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "alpha-justice",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# load thalamic hierarchy and save dataframe\n",
    "\n",
    "THhierfilename = 'hier_th_extracted_CC_TC_CT_norm.csv'\n",
    "\n",
    "with open('lib/' + THhierfilename, 'rb') as f:\n",
    "    THhier_df = pd.read_csv(f, sep=',', names=['area', 'THhier'], header=0)\n",
    "THhier_df = THhier_df.drop([34]) # remove RT area  reticular thalamic area  need to change\n",
    "THhier_df = THhier_df.reset_index(drop=True)\n",
    "print(np.array(THhier_df['THhier']).shape)\n",
    "THhier_df = pd.DataFrame(np.array(THhier_df['THhier']), columns=['THhier'], index=np.array(THhier_df['area']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "583980b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ta = np.array(THhier_df).squeeze()\n",
    "print(ta.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d91d4fe",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d201a93c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (10, 7))\n",
    "plt.scatter(range(len(thal_acr)), ta[np.argsort(ta)])\n",
    "plt.xticks(range(len(thal_acr)), np.array(thal_acr)[np.argsort(ta)], rotation = 90, size = 15)\n",
    "plt.ylabel('Thalamic hierarchy index')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f122ca99",
   "metadata": {},
   "outputs": [],
   "source": [
    "list(np.array(thal_acr)[np.argsort(ta)]).index(\n",
    "    'AD')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d533d310",
   "metadata": {},
   "outputs": [],
   "source": [
    "THhier_df.to_csv('results/dfTHhier.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed0aace4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# design using thalamocortical conections to generate a thalamic hierarchy\n",
    "plt.imshow(sortedconnthcx)\n",
    "sortedconnthcx.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a22c9a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "THhierRawfilename = 'hier_th_extracted_CC_TC_CT_norm.csv'\n",
    "\n",
    "with open('lib/' + THhierfilename, 'rb') as f:\n",
    "    THhierRaw_df = pd.read_csv(f, sep=',', names=['area', 'THhier'], header=0)\n",
    "THhierRaw_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "surprised-wilson",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6.7 ('base')",
   "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.6.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "aceb219ff324d9a10c53094f909dc059ecbb47474a96ca91c9b06a1e6187b787"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}