loopMeasurement.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "representative-associate",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:37:29.661327Z",
     "start_time": "2022-06-30T01:37:29.658260Z"
    }
   },
   "outputs": [],
   "source": [
    "# loop measurements and firing rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "weekly-albany",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:37:33.470723Z",
     "start_time": "2022-06-30T01:37:29.665094Z"
    }
   },
   "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, plotXYcomprCorrNotNormalize, plotXYcomprResidual, plotXYLogReg\n",
    "import json\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\n",
    "from collections import defaultdict\n",
    "from matplotlib.colors import TwoSlopeNorm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "simplified-portal",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:37:33.794561Z",
     "start_time": "2022-06-30T01:37:33.472618Z"
    }
   },
   "outputs": [],
   "source": [
    "# load arealist\n",
    "area_list = list(np.loadtxt('results/area_list.csv', delimiter=',', dtype='str'))\n",
    "# load connectivity \n",
    "conn_cxcx = np.loadtxt('results/conn_cxcx.csv', delimiter=',')\n",
    "# load interneuron density\n",
    "normPVgrad_df = pd.read_csv('results/dfPV.csv', header = 0, index_col=0)\n",
    "PVgrad_df = pd.read_csv('results/dfPVraw.csv', header = 0, index_col= 0)\n",
    "p = np.array(PVgrad_df)\n",
    "pvList = p/np.max(p)\n",
    "# load model\n",
    "hierarchy_df = pd.read_csv('results/dfHier.csv', header = 0, index_col=0)\n",
    "with open('results/model_V1.pkl', 'rb') as f:\n",
    "    model_V1 = pickle.load(f)\n",
    "# \n",
    "parameters = model_V1.parameters\n",
    "\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42\n",
    "plt.rcParams['svg.fonttype'] = 'none'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "together-tiffany",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:37:33.806673Z",
     "start_time": "2022-06-30T01:37:33.804472Z"
    }
   },
   "outputs": [],
   "source": [
    "noLabelAreas = ['SSp-ll','SSp-ul','SSp-n','SSp-m','SSp-tr','VISli','VISpl','VISrl','RSPd','AUDv']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hybrid-maintenance",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:37:33.815514Z",
     "start_time": "2022-06-30T01:37:33.812613Z"
    }
   },
   "outputs": [],
   "source": [
    "# ## get the connection strength for each area and test if it matches with firing rate pattern.\n",
    "# # generate normalized connection matrix\n",
    "# conn_cxcx_norm = conn_cxcx / np.max(conn_cxcx)\n",
    "# threshold = 0.2 # binary threshold\n",
    "\n",
    "# use the scaled connection matrix\n",
    "conn_sc = np.power(conn_cxcx, parameters['kSquishFln'])\n",
    "# conn_sc = conn_cxcx\n",
    "# conn_sc = conn_sc * (0.5 - 1 + pref_matrix)\n",
    "# conn_sc = conn_sc * (conn_sc > 0)\n",
    "conn_cxcx_norm = conn_sc / np.max(conn_sc)\n",
    "threshold = 0.54 # 0.2   0.6  0.56\n",
    "\n",
    "# thresholding the connection matrix\n",
    "conn_cxcx_threshold =  conn_cxcx_norm * (conn_cxcx_norm > threshold)# threholding "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "legislative-posting",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:37:33.957736Z",
     "start_time": "2022-06-30T01:37:33.817907Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(conn_cxcx_norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "responsible-example",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:49:12.785718Z",
     "start_time": "2022-06-30T01:49:12.780341Z"
    }
   },
   "outputs": [],
   "source": [
    "def makePlot(value, Xlabel, Ylabel, area_list, figureSize, fontSize, saveFig, fileName):\n",
    "    fig = plt.figure(figsize=figureSize, dpi=300, facecolor='w', edgecolor='k')\n",
    "    plt.rcParams.update({'font.size': fontSize})\n",
    "        \n",
    "    plt.plot(value)\n",
    "    plt.xticks(range(len(area_list)), area_list, rotation=90, fontsize=fontSize*0.4)\n",
    "    plt.xlabel(Xlabel)\n",
    "    # plt.grid('on')\n",
    "    plt.ylabel(Ylabel)\n",
    "\n",
    "    ax = plt.gca()\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\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": "fallen-abuse",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-07-01T00:21:34.041360Z",
     "start_time": "2022-07-01T00:21:33.272801Z"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "positive-romance",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:39:18.550338Z",
     "start_time": "2022-06-30T01:39:17.943415Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "getConnHist(conn_cxcx_norm, 'linear', threshold, (5,3), 12, True, 'loopMeasurement/plainConnection.svg')\n",
    "# getConnHist(conn_cxcx_threshold, 'log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mounted-monaco",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:39:24.285771Z",
     "start_time": "2022-06-30T01:39:23.152123Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cycDict = getNnodeLoop(conn_cxcx_norm, True, area_list, (5,3), 12, False, \n",
    "                       ['loopMeasurement/2nodeLoop.svg', 'loopMeasurement/3nodeLoop.svg'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "breathing-ladder",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:39:24.291354Z",
     "start_time": "2022-06-30T01:39:24.287926Z"
    }
   },
   "outputs": [],
   "source": [
    "cycDict_threshold = getNnodeLoop(conn_cxcx_threshold, False, area_list, (5,3), 12, False, \n",
    "                       ['loopMeasurement/2nodeLoopThres.svg', 'loopMeasurement/3nodeLoopThres.svg'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "industrial-tamil",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:39:25.727756Z",
     "start_time": "2022-06-30T01:39:25.440278Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# show the binarized connection matrix based on normalized connections\n",
    "\n",
    "binaryConn = conn_cxcx_norm > threshold\n",
    "plt.imshow(binaryConn)\n",
    "plt.figure()\n",
    "plt.imshow(conn_cxcx_norm * binaryConn)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "demographic-belly",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:39:27.396614Z",
     "start_time": "2022-06-30T01:39:27.019118Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# \n",
    "\n",
    "\n",
    "G = nx.DiGraph(binaryConn)\n",
    "# nx.draw_random(G)\n",
    "\n",
    "def getClusteringCoefficient(G, area_list):\n",
    "    ccDict = nx.clustering(G)\n",
    "    n = len(area_list)\n",
    "    ccList = []\n",
    "    for i in range(n):\n",
    "        ccList.append(ccDict[i])\n",
    "    return ccList\n",
    "cc = getClusteringCoefficient(G, area_list)\n",
    "\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(cc)\n",
    "plt.xticks(range(len(area_list)), area_list, rotation=90)\n",
    "plt.grid()\n",
    "print(cc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3f1c5d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "inputStr = np.sum(conn_cxcx_norm, axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "balanced-packaging",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:40:06.355445Z",
     "start_time": "2022-06-30T01:39:46.880133Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# compare input strength and FR\n",
    "\n",
    "X = np.sum(conn_cxcx_norm, axis = 1)\n",
    "inputStr = np.sum(conn_cxcx_norm, axis = 1)\n",
    "X = (X - X.min()) / (X.max() - X.min())\n",
    "Y = np.array(model_V1.persistentact_df['persistentact']) \n",
    "Xlabel = 'Input strength'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_FRcompr1.svg')\n",
    "inpCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_FRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/eig_FRcompr3.svg')\n",
    "inpLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state', True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_FRLogReg.svg') # Y is treated as binary with a threshold of frThreshold\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "better-fishing",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:45:29.400558Z",
     "start_time": "2022-06-30T01:45:12.915434Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# compare eigcent and FR\n",
    "eigV, eigCent = np.linalg.eig(conn_cxcx_norm)\n",
    "ind = np.argmax(eigV)\n",
    "print(ind)\n",
    "X = np.abs(eigCent[:,ind])\n",
    "eigPlain = X\n",
    "X = (X - X.min()) / (X.max() - X.min())\n",
    "Y = np.array(model_V1.persistentact_df['persistentact']) \n",
    "Xlabel = 'Eigencentrality'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/eig_FRcompr1.svg')\n",
    "eigCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/eig_FRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/eig_FRcompr3.svg')\n",
    "eigLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state',  True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/eig_FRLogReg.svg') # Y is treated as binary with a threshold of frThreshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "demanding-atlanta",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:52:40.794694Z",
     "start_time": "2022-06-30T01:52:40.787803Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# use the hierarchy scaling   start new\n",
    "# mm11\n",
    "# # use the scaled connection matrix \n",
    "# conn_sc = np.power(conn_cxcx, parameters['kSquishFln']) \n",
    "# newPref = model_V1.pref_matrix\n",
    "# conn_sc_Hiermod = conn_sc * (2*newPref - 1)\n",
    "# # conn_sc_Hiermod = conn_sc * newPref / (1 - newPref)\n",
    "# print(np.min(conn_sc_Hiermod))\n",
    "# conn_sc_Hiermod = conn_sc_Hiermod - np.min(conn_sc_Hiermod)\n",
    "\n",
    "# np.fill_diagonal(conn_sc_Hiermod, 0)\n",
    "# conn_sc_Hiermod = conn_sc_Hiermod * (conn_sc_Hiermod > 0)\n",
    "# conn_cxcx_norm_Hiermod = conn_sc_Hiermod / np.max(np.abs(conn_sc_Hiermod))\n",
    "\n",
    "# # use the scaled connection matrix\n",
    "conn_sc = np.power(conn_cxcx, parameters['kSquishFln'])\n",
    "# conn_cxcx_norm = conn_sc / np.max(conn_sc)\n",
    "newPref = model_V1.pref_matrix\n",
    "pvMat = np.matlib.repmat( np.array(normPVgrad_df), 1, len(area_list))\n",
    "# pvMat = np.matlib.repmat( pvList, 1, len(area_list))\n",
    "# hierMat = np.matlib.repmat( np.array(hierarchy_df), 1, len(area_list))\n",
    "\n",
    "# coeff = newPref - (1 - newPref)\n",
    "coeff = newPref - pvMat * (1 - newPref)\n",
    "# coeff = newPref * (1 + 1 * pvMat) - pvMat * (1 - newPref)\n",
    "\n",
    "# coeff = 1 - pvMat\n",
    "# coeff = newPref / pvMat\n",
    "# coeff = newPref / pvMat + newPref - (1 - newPref)\n",
    "# coeff = hierMat\n",
    "\n",
    "# coeff[coeff < 0] = 0\n",
    "\n",
    "minCoeff = np.min(coeff)\n",
    "print(minCoeff)\n",
    "\n",
    "# conn_cxcx_norm_Hiermod = conn_cxcx_norm * (2 * (newPref - minPref) )\n",
    "\n",
    "# conn_cxcx_norm_Hiermod = conn_cxcx_norm * (coeff - minCoeff + 0.1)\n",
    "conn_cxcx_norm_Hiermod = conn_cxcx_norm * coeff\n",
    "\n",
    "# conn_cxcx_norm_Hiermod = conn_cxcx_norm * newPref / (1 - newPref)\n",
    "print(np.min(conn_cxcx_norm_Hiermod))\n",
    "# conn_cxcx_norm_Hiermod = conn_cxcx_norm_Hiermod - np.min(conn_cxcx_norm_Hiermod)\n",
    "\n",
    "np.fill_diagonal(conn_cxcx_norm_Hiermod, 0)\n",
    "# conn_cxcx_norm_Hiermod = conn_cxcx_norm_Hiermod * (conn_cxcx_norm_Hiermod > 0)\n",
    "\n",
    "conn_cxcx_norm_Hiermod = conn_cxcx_norm_Hiermod / np.max(np.abs(conn_cxcx_norm_Hiermod))\n",
    "\n",
    "threshold = 0.2 # 0.87 #0.08  0.45 0.34\n",
    "#0.43   \n",
    "\n",
    "conn_cxcx_Hiermod_threshold =  conn_cxcx_norm_Hiermod * (conn_cxcx_norm_Hiermod > threshold) # threholding "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd7a328a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84f57d6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('results/conn_cxcx_norm.json', 'w') as f:\n",
    "    json.dump(conn_cxcx_norm.tolist(), f)\n",
    "with open('results/conn_cxcx_norm_Hiermod.json', 'w') as f:\n",
    "    json.dump(conn_cxcx_norm_Hiermod.tolist(), f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "equipped-darwin",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:52:41.612983Z",
     "start_time": "2022-06-30T01:52:41.606989Z"
    }
   },
   "outputs": [],
   "source": [
    "def makeHeatPlot(mat, area_list, 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",
    "    cbar = plt.colorbar()\n",
    "    cbar.set_label(cbarlabel)\n",
    "    plt.xticks(range(len(area_list)), area_list, rotation=90, fontsize=fontSize*0.4)\n",
    "    plt.yticks(range(len(area_list)), area_list, 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": "danish-plenty",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:52:42.827481Z",
     "start_time": "2022-06-30T01:52:42.114009Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "makeHeatPlot(coeff, area_list, 'Source area', 'Target area', 'Value','PRGn', (5, 3), 12, True, 'loopMeasurement/heat_coeff.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "treated-crazy",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:52:44.945909Z",
     "start_time": "2022-06-30T01:52:44.213682Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "makeHeatPlot(conn_cxcx_norm_Hiermod, area_list, 'Source area', 'Target area', 'Value','PRGn', (5, 3), 12, True, 'loopMeasurement/heat_modConn.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "processed-smile",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:52:45.545113Z",
     "start_time": "2022-06-30T01:52:45.231091Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "getConnHist(conn_cxcx_norm_Hiermod, 'linear', threshold, (5,3), 12, True, 'loopMeasurement/hierModConnectionHist.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "hidden-section",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:52:48.976246Z",
     "start_time": "2022-06-30T01:52:47.507957Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# cycDict_Hiermod = getNnodeLoop(conn_cxcx_norm_Hiermod, True, area_list, (10,6), 12, False, \n",
    "#                           ['loopMeasurement/2NodeLoopPVmodconn.svg', 'loopMeasurement/3NodeLoopPVmodconn.svg'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "floating-abortion",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:46:19.951320Z",
     "start_time": "2022-06-30T01:46:19.662728Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# show the binarized connection matrix based on normalized connections\n",
    "binaryConn = conn_cxcx_norm_Hiermod > threshold\n",
    "plt.imshow(binaryConn)\n",
    "plt.figure()\n",
    "plt.imshow(conn_cxcx_norm_Hiermod * binaryConn)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "retired-settlement",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-30T01:50:09.214142Z",
     "start_time": "2022-06-30T01:50:09.042568Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def makeScatter(X, Y, area_list, figTitle, xlabel, ylabel, 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.scatter(X, Y)\n",
    "\n",
    "    plt.title(figTitle)\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.ylabel(ylabel)\n",
    "    axes = plt.gca()\n",
    "    axes.spines['top'].set_visible(False)\n",
    "    axes.spines['right'].set_visible(False)\n",
    "    if saveFig:\n",
    "        fig.savefig('figure/' + fileName , dpi=300, bbox_inches='tight', transparent=True)\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "traditional-stability",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-27T05:14:08.641509Z",
     "start_time": "2022-04-27T05:13:44.760615Z"
    },
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# key test  PV vs. FR pattern\n",
    "X = np.array(normPVgrad_df).reshape((-1))\n",
    "Y = np.array(model_V1.persistentact_df['persistentact'])\n",
    "Xlabel = 'PVdensity'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/PVFRcompr1.svg')\n",
    "PVCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/PVFRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/PVFRcompr3.svg')\n",
    "PVLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state', True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/PVFRLogReg.svg') # Y is treated as binary with a threshold of frThreshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nasty-guest",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-27T05:14:26.848767Z",
     "start_time": "2022-04-27T05:14:25.624049Z"
    },
    "code_folding": [],
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# key test  Hier vs. FR pattern\n",
    "X = np.array(hierarchy_df).reshape((-1))\n",
    "Y = np.array(model_V1.persistentact_df['persistentact'])\n",
    "Xlabel = 'Hierarchy'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/HierFRcompr1.svg')\n",
    "HierCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/HierFRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/HierFRcompr3.svg')\n",
    "\n",
    "\n",
    "try:\n",
    "    HierLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state', True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/HierFRLogReg.svg') # Y is treated as binary with a threshold of frThreshold\n",
    "except Exception as e:\n",
    "    print(\"An error occurred:\", e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24dff0b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "inputStrHierMod = np.sum(conn_cxcx_norm_Hiermod, axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78d06d87",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare cell type-specific input strength and FR\n",
    "\n",
    "X = np.sum(conn_cxcx_norm_Hiermod, axis = 1)\n",
    "inputStrHierMod = np.sum(conn_cxcx_norm_Hiermod, axis = 1)\n",
    "print(X.min(), X.max())\n",
    "X = (X - X.min()) / (X.max() - X.min())\n",
    "Y = np.array(model_V1.persistentact_df['persistentact']) \n",
    "Xlabel = 'Input strength (cell type-specific)'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/inpHiermod_FRcompr1.svg')\n",
    "inpHiermodCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/inpHiermod_FRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/eig_FRcompr3.svg')\n",
    "inpHiermodLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state',  True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/inpHiermod_FRLogReg.svg') # Y is treated as binary with a threshold of frThreshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mental-possible",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-27T05:14:57.468845Z",
     "start_time": "2022-04-27T05:14:35.011258Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# key test  eigCent vs. FR pattern\n",
    "# get the eigen centrality of the conn matrix\n",
    "# mm33\n",
    "\n",
    "frThreshold = 5\n",
    "\n",
    "eigV, eigCent = np.linalg.eig(conn_cxcx_norm_Hiermod)\n",
    "# eigV, eigCent = np.linalg.eig(conn_sc)\n",
    "ind = np.argmax(eigV)\n",
    "\n",
    "X1 = np.abs(eigCent[:,ind])\n",
    "print(X1.min(), X1.max())\n",
    "eigHiermod = X1\n",
    "X = (X1 - X1.min()) / (X1.max() - X1.min())\n",
    "Y = np.array(model_V1.persistentact_df['persistentact'])\n",
    "Xlabel = 'Eigencentrality (cell type-specific)'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/eigHiermod_FRcompr1.svg')\n",
    "eigHiermodCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/eigHiermod_FRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/eigHiermod_FRcompr3.svg')\n",
    "eigHiermodLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state', True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/eigHiermod_FRLogReg.svg') # Y is treated as binary with a threshold of frThreshold\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7c3060a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newsim: compare alternative connectivity measures with FR.  sign only conn_cxcx_norm  and noPV conn_cxcx_norm.\n",
    "\n",
    "# # use the scaled connection matrix\n",
    "conn_sc = np.power(conn_cxcx, parameters['kSquishFln'])\n",
    "# conn_cxcx_norm = conn_sc / np.max(conn_sc)\n",
    "newPref = model_V1.pref_matrix\n",
    "\n",
    "# coeff = newPref - (1 - newPref)\n",
    "# change newPref matrix, all elements larger than 0.5 to 1, smaller than 0.5 to -1\n",
    "coeff_sign_only = np.ones(newPref.shape)\n",
    "coeff_sign_only[newPref < 0.5] = -1\n",
    "plt.figure()\n",
    "plt.imshow(coeff_sign_only)\n",
    "\n",
    "coeff_noPV = newPref - (1 - newPref)\n",
    "plt.figure()\n",
    "plt.imshow(coeff_noPV)\n",
    "\n",
    "# sign only conn_cxcx_norm\n",
    "conn_cxcx_norm_sign_only = conn_cxcx_norm * coeff_sign_only\n",
    "np.fill_diagonal(conn_cxcx_norm_sign_only, 0)\n",
    "conn_cxcx_norm_sign_only = conn_cxcx_norm_sign_only / np.max(np.abs(conn_cxcx_norm_sign_only))\n",
    "# noPV conn_cxcx_norm\n",
    "conn_cxcx_norm_noPV = conn_cxcx_norm * coeff_noPV\n",
    "np.fill_diagonal(conn_cxcx_norm_noPV, 0)\n",
    "conn_cxcx_norm_noPV = conn_cxcx_norm_noPV / np.max(np.abs(conn_cxcx_norm_noPV))\n",
    "\n",
    "\n",
    "#use imshow to make a heatmap of conn_cxcx_norm_sign_only matrix, add colorbar\n",
    "plt.figure()\n",
    "plt.imshow(conn_cxcx_norm_sign_only)\n",
    "plt.colorbar()\n",
    "\n",
    "#use imshow to make a heatmap of conn_cxcx_norm_noPV matrix, add colorbar\n",
    "plt.figure()\n",
    "plt.imshow(conn_cxcx_norm_noPV)\n",
    "plt.colorbar()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "unknown-january",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-22T13:44:18.902872Z",
     "start_time": "2022-06-22T13:44:18.898898Z"
    }
   },
   "outputs": [],
   "source": [
    "# Newsim: compare sign only input strength and FR \n",
    "\n",
    "X = np.sum(conn_cxcx_norm_sign_only, axis = 1)\n",
    "# inputStrHierMod??? = np.sum(conn_cxcx_norm_sign_only, axis = 1)\n",
    "print(X.min(), X.max())\n",
    "X = (X - X.min()) / (X.max() - X.min())\n",
    "Y = np.array(model_V1.persistentact_df['persistentact'])\n",
    "Xlabel = 'Input strength (sign only)'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_sign_only_FRcompr1.svg')\n",
    "inpHiermodCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_sign_only_FRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/eig_FRcompr3.svg')\n",
    "\n",
    "try:\n",
    "    inpHiermodLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state',  True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_sign_only_FRLogReg.svg') # Y is treated as binary with a threshold of frThreshold\n",
    "except Exception as e:\n",
    "    print(\"An error occurred:\", e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd81259a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newsim: compare noPV input strength and FR \n",
    "\n",
    "X = np.sum(conn_cxcx_norm_noPV, axis = 1)\n",
    "# inputStrHierMod??? = np.sum(conn_cxcx_norm_noPV, axis = 1)\n",
    "print(X.min(), X.max())\n",
    "X = (X - X.min()) / (X.max() - X.min())\n",
    "Y = np.array(model_V1.persistentact_df['persistentact'])\n",
    "Xlabel = 'Input strength (noPV)'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_nopv_FRcompr1.svg')\n",
    "inpHiermodCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_nopv_FRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/eig_FRcompr3.svg')\n",
    "\n",
    "try:\n",
    "    inpHiermodLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state',  True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/inp_nopv_FRLogReg.svg') # Y is treated as binary with a threshold of frThreshold\n",
    "except Exception as e:\n",
    "    print(\"An error occurred:\", e)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "quantitative-significance",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-07-01T18:36:05.283213Z",
     "start_time": "2022-07-01T18:36:05.278865Z"
    }
   },
   "outputs": [],
   "source": [
    "# save eigcent and input strength\n",
    "\n",
    "with open('results/inputStr.json', 'w') as f:\n",
    "    json.dump(list(inputStr), f)\n",
    "\n",
    "with open('results/inputStrHiermod.json', 'w') as f:\n",
    "    json.dump(list(inputStrHierMod), f)\n",
    "\n",
    "with open('results/hierEigCent.json', 'w') as f:\n",
    "    json.dump(list(eigHiermod), f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "electronic-victorian",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-22T13:44:02.570132Z",
     "start_time": "2022-06-22T13:44:02.563071Z"
    }
   },
   "outputs": [],
   "source": [
    "# bar plot comparing the R^2 for different settings.\n",
    "def plotBar(metric, corrList, labelList, figureSize, fontSize, saveFig, fileName):\n",
    "    fig = plt.figure(figsize=figureSize, dpi=300, facecolor='w', edgecolor='k')\n",
    "    plt.rcParams.update({'font.size': fontSize})\n",
    "    \n",
    "    yData = [c[metric] for c in corrList]\n",
    "    xData = range(len(yData))\n",
    "    plt.bar(xData, yData)\n",
    "    plt.xticks(xData, labelList)\n",
    "    plt.axhline(0, 0, 1, color='k')\n",
    "    plt.ylabel(metric)\n",
    "    axes = plt.gca()\n",
    "    axes.spines['top'].set_visible(False)\n",
    "    axes.spines['right'].set_visible(False)\n",
    "    if saveFig:\n",
    "        fig.savefig('figure/' + fileName, dpi=300, bbox_inches='tight', transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ffe3a64",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for main figure\n",
    "corrList = [inpCorr, inpHiermodCorr]\n",
    "logRegList = [inpLogReg, inpHiermodLogReg]\n",
    "labelList = ['Input', 'Input\\n (cell type\\nspecific)']\n",
    "# change font size since this is used for inset\n",
    "plotBar('r', corrList, labelList, (5,3), 24, True, 'loopMeasurement/connectivityFRCorrBar.svg') \n",
    "plotBar('prediction score', logRegList, labelList, (5,3), 24, True, 'loopMeasurement/connectivityFRLogRegBar.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "perceived-hindu",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-06-22T13:44:03.207241Z",
     "start_time": "2022-06-22T13:44:03.192516Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# for supplementary figure\n",
    "corrList = [inpCorr, inpHiermodCorr, eigCorr, eigHiermodCorr]\n",
    "logRegList = [inpLogReg, inpHiermodLogReg, eigLogReg, eigHiermodLogReg]\n",
    "labelList = ['Input', 'Input\\n (cell type\\nspecific)', 'EC', 'EC\\n (cell type\\nspecific)']\n",
    "plotBar('r', corrList, labelList, (5,3), 12, True, 'loopMeasurement/connectivityFRCorrBarSupp.svg')\n",
    "plotBar('prediction score', logRegList, labelList, (5,3), 12, True, 'loopMeasurement/connectivityFRLogRegBarSupp.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "broad-professional",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "sLoop = np.diag(conn_cxcx_norm @ conn_cxcx_norm)\n",
    "makePlot(sLoop, 'Area', 'Loop measurements', area_list, (5,3), 12, True, 'loopMeasurement/plainLoop.svg')\n",
    "\n",
    "ssLoop = np.diag(conn_cxcx_norm @ conn_cxcx_norm @ conn_cxcx_norm)\n",
    "makePlot(ssLoop, 'Area', 'Loop measurements(L3)', area_list, (5,3), 12, True, 'loopMeasurement/plainLoop3.svg')\n",
    "\n",
    "\n",
    "sLoopHier = np.diag(conn_cxcx_norm_Hiermod @ conn_cxcx_norm_Hiermod)\n",
    "makePlot(sLoopHier, 'Area', 'Loop measurements', area_list, (5,3), 12, True, 'loopMeasurement/hierLoop.svg')\n",
    "\n",
    "ssLoopHier = np.diag(conn_cxcx_norm_Hiermod @ conn_cxcx_norm_Hiermod @ conn_cxcx_norm_Hiermod)\n",
    "makePlot(ssLoopHier, 'Area', 'Loop measurements(L3)', area_list, (5,3), 12, True, 'loopMeasurement/hierLoop3.svg')\n",
    "\n",
    "# compare length 2 and length 3 loop\n",
    "makeScatter(sLoopHier, ssLoopHier, area_list, 'Length 2 loop vs Length 3 loop', 'Length 2 loop', 'Length 3 loop', \n",
    "    (5, 3), 12, True, 'loopMeasurement/l2vsl3.svg')\n",
    "r, pvalue = sp.stats.pearsonr(sLoopHier, ssLoopHier)\n",
    "print('correreation between length 2 and length 3 loop: ', r, pvalue)\n",
    "# save the loop \n",
    "with open('results/plainLoopStr.json', 'w') as f:\n",
    "    json.dump(list(sLoop), f)\n",
    "\n",
    "with open('results/plainLoopStr3.json', 'w') as f:\n",
    "    json.dump(list(ssLoop), f)\n",
    "\n",
    "with open('results/hierLoopStr.json', 'w') as f:\n",
    "    json.dump(list(sLoopHier), f)\n",
    "\n",
    "with open('results/hierLoopStr3.json', 'w') as f:\n",
    "    json.dump(list(ssLoopHier), f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "410a2506",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare raw loop strength and FR\n",
    "# the result is not good. so we start to consider the loop strength to predict core areas.\n",
    "\n",
    "X = sLoop\n",
    "print(X.min(), X.max())\n",
    "X = (X - X.min()) / (X.max() - X.min())\n",
    "Y = np.array(model_V1.persistentact_df['persistentact']) \n",
    "Xlabel = 'Loop strength'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/loop_FRcompr1.svg')\n",
    "loopHierCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/loop_FRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/eig_FRcompr3.svg')\n",
    "loopHierLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state',  True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/loop_FRLogReg.svg') # Y is treated as binary with a threshold of frThreshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5438c73c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compare cell type-specific loop strength and FR\n",
    "# the result is not good. so we start to consider the loop strength to predict core areas.\n",
    "\n",
    "X = sLoopHier\n",
    "print(X.min(), X.max())\n",
    "X = (X - X.min()) / (X.max() - X.min())\n",
    "Y = np.array(model_V1.persistentact_df['persistentact']) \n",
    "Xlabel = 'Loop strength (cell type)'\n",
    "Ylabel = 'Firing rate (Hz)'\n",
    "frThreshold = 5\n",
    "plotXYcompr(X, Y, Xlabel, Ylabel, area_list, (5, 3), 12, True,\n",
    "            'loopMeasurement/loopHier_FRcompr1.svg')\n",
    "loopHierCorr = plotXYcomprCorrNotNormalize(X, Y, Xlabel, Ylabel, parameters, False, area_list, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/loopHier_FRcompr2.svg')\n",
    "# plotXYcomprResidual(X, Y, Xlabel, frThreshold, (5, 3), 12, True,\n",
    "#             'loopMeasurement/eig_FRcompr3.svg')\n",
    "loopHierLogReg = plotXYLogReg(X, Y, Xlabel, 'Delay period state',  True, parameters, area_list, noLabelAreas, frThreshold, (5, 3), 12, True,\n",
    "            'loopMeasurement/loopHier_FRLogReg.svg') # Y is treated as binary with a threshold of frThreshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e57633af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newsim: conn_cxcx_norm_sign_only    conn_cxcx_norm_noPV\n",
    "\n",
    "sLoop_sign_only = np.diag(conn_cxcx_norm_sign_only @ conn_cxcx_norm_sign_only)\n",
    "makePlot(sLoop_sign_only, 'Area', 'Loop measurements', area_list, (5,3), 12, True, 'loopMeasurement/sign_onlyLoop.svg')\n",
    "\n",
    "ssLoop_sign_only = np.diag(conn_cxcx_norm_sign_only @ conn_cxcx_norm_sign_only @ conn_cxcx_norm_sign_only)\n",
    "makePlot(ssLoop_sign_only, 'Area', 'Loop measurements(L3)', area_list, (5,3), 12, True, 'loopMeasurement/sign_onlyLoop3.svg')\n",
    "\n",
    "sLoop_noPV = np.diag(conn_cxcx_norm_noPV @ conn_cxcx_norm_noPV)\n",
    "makePlot(sLoop_noPV, 'Area', 'Loop measurements', area_list, (5,3), 12, True, 'loopMeasurement/noPVLoop.svg')\n",
    "\n",
    "ssLoop_noPV = np.diag(conn_cxcx_norm_noPV @ conn_cxcx_norm_noPV @ conn_cxcx_norm_noPV)\n",
    "makePlot(ssLoop_noPV, 'Area', 'Loop measurements(L3)', area_list, (5,3), 12, True, 'loopMeasurement/noPVLoop3.svg')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "002bbc36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save the loop \n",
    "\n",
    "with open('results/sign_onlyLoopStr.json', 'w') as f:\n",
    "    json.dump(list(sLoop_sign_only), f)\n",
    "\n",
    "with open('results/sign_onlyLoopStr3.json', 'w') as f:\n",
    "    json.dump(list(ssLoop_sign_only), f)\n",
    "\n",
    "with open('results/noPVLoopStr.json', 'w') as f:\n",
    "    json.dump(list(sLoop_noPV), f)\n",
    "\n",
    "with open('results/noPVLoopStr3.json', 'w') as f:\n",
    "    json.dump(list(ssLoop_noPV), f)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c50f16d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "newList = [area_list.index(x) for x in ['FRP','VISa','VISam','MOs','ACAd']]\n",
    "print(newList)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b14c8f0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "partOfConn = conn_cxcx_norm_Hiermod[newList, :][:, newList]\n",
    "print(partOfConn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8f0df2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sLoopHierPart = np.diag(partOfConn @ partOfConn)\n",
    "plt.plot(sLoopHierPart)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a8895ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeHist(X, area_list, nBins, logScale, figTitle, xlabel, ylabel, figureSize, fontSize, saveFig, fileName):\n",
    "    fig = plt.figure(figsize=figureSize, dpi=300, facecolor='w', edgecolor='k')\n",
    "    plt.rcParams.update({'font.size': fontSize})\n",
    "    \n",
    "    if logScale:\n",
    "        hist, bins = np.histogram(X, bins=nBins)\n",
    "        logbins = np.logspace(np.log10(bins[0]),np.log10(bins[-1]),len(bins))\n",
    "        plt.hist(X, bins=logbins)\n",
    "        plt.xscale('log')\n",
    "    else:\n",
    "        plt.hist(X, bins=nBins, color='#0504aa', alpha=0.7, rwidth=0.85)\n",
    "\n",
    "    plt.title(figTitle)\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.ylabel(ylabel)\n",
    "    axes = plt.gca()\n",
    "    axes.spines['top'].set_visible(False)\n",
    "    axes.spines['right'].set_visible(False)\n",
    "    if saveFig:\n",
    "        fig.savefig('figure/' + fileName , dpi=300, bbox_inches='tight', transparent=True)\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "171e78d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # histogram of loop 2 strength\n",
    "loopStrList = []\n",
    "for i in range(len(area_list)):\n",
    "    for j in range(len(area_list)):\n",
    "        if i==j:\n",
    "            continue\n",
    "        loop = conn_cxcx_norm[i, j] * conn_cxcx_norm[j, i]\n",
    "        loopStrList.append(loop)\n",
    "makeHist(loopStrList, area_list, 100, False, 'Distribution of Length 2 Loops', 'Loop measures', 'Frequency',\n",
    "    (5, 3), 12, True, 'loopMeasurement/Rawl2Hist.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d972192",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # histogram of loop 3 strength\n",
    "loop3StrList = []\n",
    "for i in range(len(area_list)):\n",
    "    for j in range(len(area_list)):\n",
    "        if i==j:\n",
    "            continue\n",
    "        for k in range(len(area_list)):\n",
    "            if k == i or k == j:\n",
    "                continue\n",
    "            loop = conn_cxcx_norm[i, j] * conn_cxcx_norm[j, k] * conn_cxcx_norm[k, i]\n",
    "            loop3StrList.append(loop)\n",
    "makeHist(loop3StrList, area_list, 100, False, 'Distribution of Length 3 Loops', 'Loop measures', 'Frequency',\n",
    "    (5, 3), 12, True, 'loopMeasurement/Rawl3Hist.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75e506db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# histogram of loop 4 strength\n",
    "loop4StrList = []\n",
    "for i in range(len(area_list)):\n",
    "    for j in range(len(area_list)):\n",
    "        if i==j:\n",
    "            continue\n",
    "        for k in range(len(area_list)):\n",
    "            if k == i or k == j:\n",
    "                continue\n",
    "            for l in range(len(area_list)):\n",
    "                if l == i or l == j or l == k:\n",
    "                    continue\n",
    "                loop = conn_cxcx_norm[i, j] * conn_cxcx_norm[j, k] * conn_cxcx_norm[k, i]\n",
    "                loop4StrList.append(loop)\n",
    "makeHist(loop4StrList, area_list, 100, False, 'Distribution of Length 4 Loops', 'Loop measures', 'Frequency',\n",
    "    (5, 3), 12, True, 'loopMeasurement/Rawl4Hist.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b36df11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# histogram of Hierarchy modified loop 2 strength\n",
    "HiermodloopStrList = []\n",
    "for i in range(len(area_list)):\n",
    "    for j in range(len(area_list)):\n",
    "        if i==j:\n",
    "            continue\n",
    "        loop = conn_cxcx_norm_Hiermod[i, j] * conn_cxcx_norm_Hiermod[j, i]\n",
    "        HiermodloopStrList.append(loop)\n",
    "makeHist(HiermodloopStrList, area_list, 100, False, 'Distribution of Cell Type-Specific Length 2 Loops', 'Loop measures', 'Frequency',\n",
    "    (5, 3), 12, True, 'loopMeasurement/hiermodl2Hist.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "736abc8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # histogram of Hierarchy modified loop 3 strength\n",
    "Hiermodloop3StrList = []\n",
    "for i in range(len(area_list)):\n",
    "    for j in range(len(area_list)):\n",
    "        if i==j:\n",
    "            continue\n",
    "        for k in range(len(area_list)):\n",
    "            if k == i or k == j:\n",
    "                continue\n",
    "            loop = conn_cxcx_norm_Hiermod[i, j] * conn_cxcx_norm_Hiermod[j, k] * conn_cxcx_norm_Hiermod[k, i]\n",
    "            Hiermodloop3StrList.append(loop)\n",
    "makeHist(Hiermodloop3StrList, area_list, 100, False, 'Distribution of Cell Type-Specific Length 3 Loops', 'Loop measures', 'Frequency',\n",
    "    (5, 3), 12, True, 'loopMeasurement/hiermodl3Hist.svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "voluntary-niagara",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get which areas show strongest loopstrength\n",
    "loopSortedAreaList = [area_list[i] for i in np.argsort(sLoop)]\n",
    "print(loopSortedAreaList[::-1])\n",
    "# get which areas show strongest loopstrength\n",
    "loopSortedAreaList = [area_list[i] for i in np.argsort(sLoopHier)]\n",
    "print(loopSortedAreaList[::-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "leading-clerk",
   "metadata": {},
   "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
}