RBC_subtypes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of possible RBC subtypes\n",
    "The code in this notebook analyses differences between the RBCs contacting cones and rods and those contacting only rods in order to find evidence for or against two types of RBCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import itemfreq\n",
    "from scipy.io import loadmat\n",
    "from scipy.spatial import ConvexHull\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from PIL import Image\n",
    "from PIL import ImageDraw\n",
    "from sklearn.mixture import GMM\n",
    "from shapely.geometry import Polygon\n",
    "%matplotlib inline\n",
    "matplotlib.rcParams.update({'font.size': 14})\n",
    "matplotlib.rcParams['font.family'] = 'sans-serif'\n",
    "matplotlib.rcParams['font.sans-serif'] = ['Arial']\n",
    "sns.set_style(\"whitegrid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cone_contacts=pd.read_pickle('data/cone_contact_predictions')\n",
    "rod_contacts=pd.read_pickle('data/rod_contact_predictions')\n",
    "BC_IDs=np.loadtxt('data/BC_IDs_new')\n",
    "ID_table=np.loadtxt('data/cell_IDs.csv',dtype=int,delimiter=',')\n",
    "BC_in_rod_area=np.loadtxt('data/BC_in_rod_area')\n",
    "BC_excluded=np.array([691,709,827,836])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rbc=BC_IDs[(BC_IDs[:,4]>70)&np.in1d(BC_IDs[:,0],BC_in_rod_area)&np.in1d(BC_IDs[:,0],BC_excluded,invert=True),0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rbc_rod_contacts=rod_contacts.ix[(rod_contacts['prediction']==1)]\n",
    "rbc_rod_contacts=rbc_rod_contacts[np.in1d(rbc_rod_contacts['cell'],rbc)].reset_index().drop('index',axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rbc_cone_contacts=cone_contacts.ix[(cone_contacts['prediction']==1)]\n",
    "rbc_cone_input=np.unique(rbc_cone_contacts[np.in1d(rbc_cone_contacts['cell'],BC_excluded,invert=True)&\\\n",
    "                                           np.in1d(rbc_cone_contacts['cell'],BC_IDs[BC_IDs[:,4]==71,0])]['cell']).astype(int)\n",
    "rbc_rod_input=BC_IDs[(BC_IDs[:,4]==71)&np.in1d(BC_IDs[:,0],rbc_cone_input,invert=True)&np.in1d(BC_IDs[:,0],BC_excluded,invert=True),0].astype(int)\n",
    "rbc_cone_contacts=rbc_cone_contacts[np.in1d(rbc_cone_contacts['cell'],rbc)].reset_index().drop('index',axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#mean(on_sac)=6327, mean(off_sac)=5403\n",
    "x_start=-5403*(0.62-0.28)/(6327-5403)+0.28\n",
    "x_end=10000*(0.62-0.28)/(6327-5403)+x_start\n",
    "x_scale=np.linspace(x_start*100,x_end*100,10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "skeleton_data=loadmat('data/skeletons_OPL_final.mat')\n",
    "skeleton_ids=skeleton_data['kn_allSkeletons_e2006_IDs'].flatten()\n",
    "skeletons=skeleton_data['kn_allSkeletons_e2006'].flatten()\n",
    "M=np.loadtxt('data/coordinate_rotation')\n",
    "soma_data=loadmat('data/soma_positions.mat')\n",
    "soma_pos=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_allSomata']\n",
    "soma_internal_ids=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_cellIDs']\n",
    "soma_line_ids=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_cellIDs_pure_forSomata']\n",
    "soma_pos=np.dot(M,soma_pos[:,:3].T).T\n",
    "soma_positions=[]\n",
    "for i in range(BC_IDs.shape[0]):\n",
    "    soma_positions.append(soma_pos[soma_line_ids[0,np.where(soma_internal_ids==BC_IDs[i,1])[1][0]]-1,:])\n",
    "soma_positions=np.array(soma_positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Number of rod contacts (Supp. figure 3B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stat_bc_contacts=pd.DataFrame(rbc,columns=['cell'])\n",
    "contact_freq_rbc=itemfreq(rbc_rod_contacts['cell'].as_matrix())\n",
    "for i in range(stat_bc_contacts.shape[0]):\n",
    "    stat_bc_contacts.loc[i,'count']=0\n",
    "    try:\n",
    "        stat_bc_contacts.ix[i,'count']=contact_freq_rbc[contact_freq_rbc[:,0]==stat_bc_contacts.ix[i,'cell'],1]\n",
    "    except ValueError:\n",
    "        continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in range(stat_bc_contacts.shape[0]):\n",
    "    if stat_bc_contacts.ix[i,'cell'] in np.unique(rbc_cone_contacts['cell']):\n",
    "        stat_bc_contacts.loc[i,'cone_contacts']=1\n",
    "    else:\n",
    "        stat_bc_contacts.loc[i,'cone_contacts']=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHMAAAB4CAYAAADBhM0lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACpVJREFUeJztnX9M1GUcx193xy4670wF9MCAkGkHFtIWY822IHIm0SYr\nhq0yltpVNtY8Z8M13N0Sslxb5cJMp8w/8tBlm9X6RbWWKyfQIoOD1omIdcQhyZ3Q9Ivf/mhcusTj\ngO+X43vP669jd9/neY/Xnu99n4cPz6OTZVlGoAn00x1AMHUImRpCyNQQQqaGiGqZkiTR09ODJEnT\nHWVGENUyfT4fRUVF+Hy+6Y4yI4hqmYLIiFOy8c2bN3P//ffzxx9/cO7cOYLBIFVVVcydO1fJbmMW\nxUbmgQMHmDVrFgBNTU1UV1fzyCOP0NDQoFSXMY8iMr/66issFgu5ublcuXKFhIQEAKxWK319fUp0\nKUCh2+yxY8e45ZZb8Hq9AKER6vP5mD9//g2vdbvduN1uAC5dujSluf7++2/i4+OntM1oyqBTcm32\nww8/xGg00t/fj9frZXBwEKfTidlsHtf1PT09FBUV0djYyK233jqpLH19fdx22210dXWRlJQ0qbai\nNYOiD0CrV69WsvmICAaDDA0NEQwGp02m0hkUlSn4D4/Hg8fjAeCnn37C7/djsViw2WxT1oeQqQIe\nj4esrKzQz6WlpaHX7e3tUyZULBqoQCAQmNB7kSJkagghU0MImSpgsVgm9F6kiAcgFbDZbLS3t+Px\neCgtLeXo0aMsXLhQPM3OVGw2GzfddBMAy5YtIyMjY8r7ELdZFTGbzZhMpnGvgEWKGJkqkpSURH9/\nv2Jrs2JkqoySC/1CpoaIidusGuui0YDmZaq1LhoNaP42q9a6aDSg2Mg8c+YMb775JvPmzWPp0qWc\nP39eFHUpjGIyg8EgDoeDBQsWYLfbMRqN1NXVceLECRoaGrDb7Up1HbMoJnPp0qX09vZit9vJz8+n\nu7sbCF/UpWQNkNZRTGZbWxspKSns27ePyspKRkZGgPBFXeXl5ZSXlwP/1QBNBrUWuaMBxWRKkkR1\ndTVWq5XU1FSsVitOpzNU1KUWai1yRwOKyczJyeGtt95SqvmIUGOROxrQ/NQklhAyNYSQqSGETA0R\nkcze3l6lcgimgLBPs++88w7Jycn09PTQ1NREeno6LpdLjWyCCAk7Mr/77jtKS0tpbm6mvr6erq4u\nFWIJJkJYmbIss3//frKzs2ltbWV4eFiNXIIJEFbm9u3bMRqNVFZW0tPTw+uvv65GLsEEGPM7c9eu\nXdf8vHfvXgC8Xi8vvPCCsqkEE2LMkWmz2bDZbLS2tmIwGMjNzSU+Pp6Ojg418wkiYMyR+cADDwBw\n8OBBnnvuOQDuvfde1q5dq04yQcSEnZrMnj2bt99+G5vNxs8//0xiYqIauaYcpQuQo4GwD0BvvPEG\nt99+O11dXWRnZ8/YB6DRAuTp+hd4NQg7Mru6unC73QwMDDBr1izS0tLIzs5WI9uUM907jShNWJku\nl4uamhpSU1Pp7u5my5YtHDp06IbX/Pjjjxw6dAiz2cy8efOIj48XxVwqEFbm5cuXSU1NBSAtLQ2d\nThe20cHBQbZt24bJZGLdunWimEslwsrMz89n06ZN5OTk0NraSl5eXthG77vvPgB2797Nww8/TFNT\nEzC+HbpEQdckkMPw9ddfy21tbfLHH38snzp1KtzHZVmW5WAwKG/dulX+/vvvZUmS5I0bN8qyLMs/\n/PCD/O67746rDVmW5bNnz8pLliyRz549O+5rYpmwO3Q9+eST1NfXo9eP/69lW7dupbu7m5SUFAwG\nA9nZ2dO+Q1csEPY229fXx/Lly0lOTkan06HT6Thy5MgNr6mpqZmygILxE1bmp59+qkYOwRQgykY0\nhJCpIcaU2djYyO+//x4qERGlItHPmN+ZkiTR0NDAZ599xuDgIC0tLeTl5bFkyRIyMzPVzCgYJ2OO\nzLi4OMrKysjMzGTnzp3Mnz8fk8lEY2OjmvkEETCmTEmSOHz4MF6vl5deeok///yTQCBAQUGBivEE\nkTCmzJUrV/Liiy9y1113sWPHDhYvXozZbOabb75RMZ4gEhTdo32yiBWgyBBTEw0hZGoIIVNDCJka\nQsjUEEKmhlBU5pkzZ0J71e3btw+Xy8WWLVsYGBhQstuYRTGZfr+fI0eOYDKZuHTpEidPnhTHLiqM\nYlvHJCYm4nA42LBhAxcuXAhVwosdupRDlS1KExIS+OuvvwD1d+iKJVSRqdfryc/Pn5YdumIJsTar\nIcTUREMImRpCyNQQQqaGEDI1hJCpIYRMDSFkagghU0MImRpCyNQQQqaGEDI1hJCpIYRMDaHaYai9\nvb3s2LGDOXPmkJmZyeOPP65W1zGDajLdbjdr164lNzeXZ555hjVr1mAwGMZ9vSRJ+Hw+BRPOLKxW\nK3Fx1+pTTabf7yc5ORn4d9vTQCDAnDlz/ve5sQq6fD6fqAe6iutVX6gmMyUlBZ/Px4IFCxgcHGT2\n7NnX/dzVBV2jo9FqtQKI/9q+itHfydWoVgPk9/upra3FbDZzxx13UFZWpka3MUVUF3QJIkNMTTSE\nkKkhhEwNodrTrBIcPXqUY8eOkZaWxsjICENDQ9TW1mI0Gq/7+ZqaGioqKkhJSYmKPJFkcrvdoaf8\nsZjRMgEeffRRiouLAXA6ndTX19PR0YHRaOSee+6hsLCQbdu2kZCQwMmTJ6moqGDXrl1cvHgRn8+H\n0+kcc5o02Twul4vffvuNgYEBPvjgg0ll+uKLL7Qv8/3336exsZHTp0+zatUqmpubqaurQ6fTUVFR\nwfDwMCUlJRQWFlJdXQ3AuXPnyMjIoLCwkJtvvlmxPMXFxWRlZWG329m9e/eEMn300UccP36czs5O\nqqqqePrpp1m8ePF1+57xMh977DGKi4s5fPgwHR0dXLlyJfSeTqfDYDAwMjICEFr+WrNmDXFxcRw8\neJDh4eHQnvJTnaezsxP496TDiWYqKSmhpKSEDRs2UFtbe8O+Z7zMUcrKynjllVcoLi6mqqoKs9lM\nWVkZBQUFuFwumpubOXXqFABffvklgUAAnU7HokWLFM3zySef8NRTT00603vvvRe2T7FooCHE1ERD\nxITM9evXT3cEVdDEd6bH42HPnj1YLBbS09P5/PPPKSgooK2tje3bt6PT6Th9+jQHDhzA6XRSV1fH\n8uXLycnJUSRPZ2cn+/fvR5Zl0tPT+fXXX0PZVq5cicPhoLCwkF9++YWamhq+/fZbWlpaGBoaYsWK\nFcydO5eGhgZkWaaoqCh0/GU4NDEy9+zZg9PpxOl00tLSQlJSEs8++yxZWVmhw1szMjIIBAKcP3+e\n9vZ2xUSO5nn55Zd59dVXaWpquibbxYsXSU9Px263k5WVhcfjYe/evZhMJhITEzl+/Dh+v5/Lly9T\nVFTEnXfeOe5+NSHzanQ6HRaLBQCDwXDNVKW8vByHw8GDDz6oaAZJkkKvT5w4EXqt1+uRZTk0jxzN\np9fr2bRpExs3bmTZsmVkZmby/PPPc+HCBV577bVx96uJ2+z69etxuVwkJCSQl5eH1+sF+N8hdPn5\n+ezcuVNxmevWrcPlcqHX66msrAxlu/vuu685bWn00J8nnngCh8MB/Dul6e/vp76+noULF47r7LVQ\ne7EyNZEkic2bN7NixQoeeuih6Y6jCDEjMxbQ3HdmLCNkagghU0MImRpCyNQQ/wBcX9m8VZxOuQAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb6e066ad68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set(font='Arial',style='white',context='paper',rc={\"xtick.major.size\": 0, \"ytick.major.size\": 4})\n",
    "with matplotlib.rc_context({\"lines.linewidth\": 0.7}):\n",
    "    plt.figure(figsize=(3/2.54,4/2.54))\n",
    "    ax=sns.pointplot(x='cone_contacts',y='count',data=stat_bc_contacts,ci=95,order=[0,1], \\\n",
    "                     linestyles='',markers='s',scale=1.5,color='black')\n",
    "    ax.set(ylabel='# rods',ylim=(0,42),yticks=[0,10,20,30,40],xlabel='',xticklabels=['Rods\\nonly','Rods +\\ncones'],xticks=[-0.1,1.1])\n",
    "    sns.despine()\n",
    "    ax.spines['left'].set_position(('outward',3))\n",
    "#     plt.savefig('figures/rbc_rod_contact_comparison.svg',bbox_inches='tight',dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### AII connectivity (Supp. figure 3C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "connectivity=loadmat('data/Helmstaedter_et_al_SUPPLinformation5.mat')['kn_allContactData_Interfaces_duplCorr_output_IDconv']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rbc_ids=ID_table[ID_table[:,3]==71,0]\n",
    "ac_ids=ID_table[ID_table[:,3]==24,0]\n",
    "rbc_ac_contacts=connectivity[np.in1d(connectivity[:,0],ac_ids)&np.in1d(connectivity[:,1],rbc_ids)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ac_connectivity=pd.DataFrame(BC_IDs[(BC_IDs[:,4]==71)&np.in1d(BC_IDs[:,0],BC_excluded,invert=True),0],columns=['rbc'])\n",
    "for i in range(ac_connectivity.shape[0]):\n",
    "    if ac_connectivity.ix[i,'rbc'] in rbc_cone_input:\n",
    "        ac_connectivity.loc[i,'cone_input']=1\n",
    "    else:\n",
    "        ac_connectivity.loc[i,'cone_input']=0\n",
    "    ac_connectivity.loc[i,'area']=np.sum(rbc_ac_contacts[rbc_ac_contacts[:,1]==ac_connectivity.ix[i,'rbc'],2])\n",
    "ac_connectivity=ac_connectivity[ac_connectivity['area']>0].reset_index().drop('index',axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHkAAAB4CAYAAADWpl3sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADh1JREFUeJztnXtM1eUfx1/fw0VFMARkCCloq8AL2IXM1AkDZxK2XCGa\npiR4CZ1bYibq0HMSdHbZyjbJS3lZS9DQTcgWupFJzRArK+WyOBCixwAvwIFC8Pv7wzg/yQvny/l+\nuZ3ntbF5bs/zdu/v5Xmez+fzfCVZlmUEfRpddwsQaI8w2Q4QJtsBwmQ7oFea3NLSwsWLF2lpaelu\nKb2CXmmyyWQiIiICk8nU3VJ6Bb3SZIEyhMl2gDDZDnC05ktTpkxhyJAh9/xMlmVqa2vJy8tTU5dA\nRawyefLkyWzatOm+n69fv96qzn766ScOHDiAq6srHh4eXLp0idbWViRJYs6cOQQHB1unWqAIq0x+\n5513ALh58yZOTk53ff6gA+BO6urq2LBhAy4uLsTHx1NXV8fYsWPR6XQ8+uijCmQLlCB1FKAwGAwA\nzJo1i7y8PJYuXWpzp+np6fj4+ODn50doaCh5eXkUFxezZMmS+/4mIyODjIwMAJqbmyktLeXEiRM8\n/PDDNuv5+++/6d+/v83t9FQ6HHj5+vqSkpLCsWPHuHr1qk2dmc1m1q1bx7hx44iMjKS0tBQAd3f3\nDhc2YmNjycrKIisri/T0dJt03El1dTWenp5UV1er1mZPo0OT3dzc+Ouvv1i2bBnXrl2zqbPU1FQq\nKirIysoiNTUVo9FIamoqmZmZvPrqqza13VkaGhpobGykoaGhW/rvEuReSGVlpfzYY4/JlZWVNrdV\nVlYmA3JZWZkKynomVg282ti/fz9Hjx7l1q1byLKMJEkcOnRIo8NPoBaKTD58+DAHDhzA2dlZKz0C\nDVBk8qRJkygsLMTf39/ynq+vr+qiBOqiyOSqqip27drFwIEDAZAkiQ8//FATYQL1UGTy9evX2b17\nt1Za7Bot5+qKTB44cCAbN25stzo1d+5c1UXZG9XV1QQEBFBeXn7fGIEtKDI5PDy83WtJklQV09UU\nFRVRVFQEwC+//MI///xDYGBgl+u4c67e7SZ7eHioLqC7KCoqIigoyPJ65syZAFy4cKFLje6KA02R\nyb/++isAra2tlJeXc/nyZaZMmaKqoK6ivr5e0fta0FUHmiKTly9f3u61uB/bRlcdaIpM/vzzz4Hb\niQLl5eWYzWbLe9YY/t94cv/+/amqqqKhoYHk5GQGDx7cif+CoCOsMvmHH35gwoQJuLi4WN4bPXo0\no0ePBrA6Nfa/8WRnZ2e2b9/O6dOnyczMfGCoUdB5rDJ548aNzJ8//76f79u3j5iYmA7babt/p6en\nM2PGDM6cOQOAj49Ph6G+/8aTbcXNzU3R+1rQVRqsMrmjRAFrEwnMZjNpaWnMmDGD0NBQjh8/DtzO\no/b29n7gb2NjY4mNjQXg4sWLREREWNXn/QgMDOTChQsUFRUxc+ZMDh8+TGBgYJeOrLtKQ4eZIWqy\ndu1a/vzzT3x9fXFwcGDUqFGUlZVRV1eHXq/H1dXVqnbaTFYjM8RoNDJy5EjKysoYMWKETW31VA2K\nBl62kpaW1pXdCf5F5F33AFxdXXFxcbH6SqYURWdyQ0MDhYWFtLS0IMsyNTU1zJ49WxNh9sSQIUOo\nra3tGQGKxMREgoOD+fnnn/Hw8MDd3V0TUfaIltmiii7Xjo6OrFq1imHDhvHRRx9x5coVrXQJVESR\nyQMGDKCwsJBbt25x8uRJKioqtNIlUBFFJm/atAlnZ2dWrFjBd999x7p167TSJVARRffkhx56iFOn\nTvHjjz8SFRVFQECARrK6Dq1Htj0BRWdycnIytbW15Obm0tjYSFJSkla6uoy2ka0WwfqegiKTTSYT\ncXFxODs7M3HiRFXWkHsCfbkOChSa7O7uzsGDB2lsbOTYsWMiNNhLUGRyamoqTU1NjBkzBpPJxJYt\nWzrVaUVFhSULYu3ataxZs4bk5GTOnTvXqfYED0bRwOuNN95g//79NnVYU1PDoUOHLLHpkpISgoOD\nRY2yhiiKQr322msEBgbi7+9vydTsbArQokWL2LlzJwUFBVbVKGtZn9zXUXQmR0ZG4ujoSGtrK7Is\n069fv053LMsyZrOZ0tJSQkNDO6xRVjuebE8oMvn48eOEhISosnYtSRIDBw6krKyM1NRUzGYzb731\nVqfbEzwAJXWur7/+uizLsrxmzRpZlmV58eLFqtXQKkHN+mR7QKxd2wGK166dnJzE2nUvQ9E9efDg\nwZYFEGFw70Gk/9gBikxuy5Nu4/vvv1dVjEAbrLpc5+Xlcfr0ab7++mumT58O3J7nnjp1iqNHj2oq\nUGA7Vpk8btw4XFxcqKqqIiwsDLg9z12wYIGW2gQqYZXJ7u7uPPPMM4SGhlJdXY23tzfFxcX4+Pho\nrU+gAoruyWvXriU3NxeA3Nxcq3fHFXQvikw2Go2WgMTy5cvFYkgvQdE8edCgQRw6dIhRo0Zx/vz5\nPp0X1ZdQdCa///77NDQ0cPDgQerq6nj33Xc71emdSQO7d+/GYDCwevVqmzdoFdwbRSZLkkRAQACT\nJk1i+PDhZGdnK+7wzqSB5uZmCgoKSElJ4eWXXyYzM1Nxe4KO6fIyGS8vL5KSkli0aBE3btzAy8sL\n6LgQXe0idHuiW8tkPD09uX79OtBxIbpWm5rbA4rOZLVDjTqdjvHjx6PX6y2F6AL1UZTjVVBQQP/+\n/fHw8GDPnj0899xzd+3S1xWoudOAPdDptWtHR0c++OCDbjFZoAzFa9fh4eGWXevF2nXvQNHa9eOP\nP05+fr5ldFtVVcVLL72kqUCB7SgaXS9btozi4mLMZrPlT9DzUTS67tevH2+++aZWWgQaochkNzc3\nsal5L0SRyW3bJkqShIKZl6CbUXRPfvHFF9HpdJSXl9Pa2kp0dLRWugQqosjkdevWUVlZyVNPPYXJ\nZGL16tVa6RKoiOJHBrXVJE+ePFncj3sJikx2cnLi1KlTjB07lnPnzqmyDcOlS5dITEwkKCgIb29v\nMXrXAEUmb9myhW3btrF3716am5s7nTRwJwUFBZZNWZ544gmb2xPcjaIARXJyMmPGjGHu3Lls27aN\nK1euWP0U9PthNBpxc3PD09OTuLg4Pv30UxwcHO76nihCtwElJZCxsbHtXs+bN8/mssqcnBy5qqpK\nlmVZTkxMlJubmzv8jShdVUa3J/L5+/uzdetWvLy8CAsLw8nJyeY2Be1RdLmur6/nyy+/pKKigmHD\nhjFr1qxuydgU8WRlKF7WjIuL00iKQCtE6aodIEy2A6y6XK9YseKuJ6zK/2aHiIdk93ysMvntt9/W\nWodAQ6wy2c/PT2sdAg3p9OW6DXG57vl06nJ98+ZN8vLyKCgo0ESUQF0UXa7PnTvHkSNHOHHiBFOn\nTmXhwoWaihOog1Umf/LJJ3z77beMHDmSadOmUV5eLnYZ6EVYNU/Oz89n6NChPPvsszz55JP3jBJ1\nlitXrrBy5UoMBoPlgdsCdbHqTN63bx+XL18mJyeH+Ph4Kisryc7OJiwszOa164yMDObPn8+4ceNY\nvHgxs2fPVnQQtbS0YDKZbNLQl/Dx8cHRsb2tVq9dDx06lISEBBISEjAajeTk5DBv3jyOHDlik6ia\nmhqGDh0K3I5y1dfX37Pu+X71ySaTSex9fQf3Ctp06fOT70V6ejoTJkwgJCSExYsXk56ejk734LtI\n29nbtsWUOJP/z73O5G43uaamhs2bN+Pq6sqYMWOIiYnpTjl9km43WaA9IgplBwiT7QBhsh2gKP2n\nt3D48GGOHj3K8OHDaW1tpbGxkc2bN+Ps7HzP76elpREXF4evr2+3a1GiJyMjw/IYpQfRJ00GeOWV\nV4iKigJAr9ezd+9eiouLcXZ2ZsKECYSHh7NhwwY8PT0pKCggLi6Ojz/+GLPZjMlkQq/XM2jQINW1\nGAwG/vjjD65du0ZWVpZNenJzc+3b5C+++IITJ05gNBqZPn06hYWFbN++HUmSiIuLo6mpiejoaMLD\nw0lJSQFu13qNGDGC8PBwBgwYoImWqKgogoKCWLJkCenp6Z3Sk52dTX5+PiUlJSQnJ7Nw4cIHPgKx\nz5o8Z84coqKiOHjwIMXFxdy6dcvymSRJODg40NraCmBZPJg9ezaOjo7s37+fpqYmSz22mlpKSkoA\n2tV3K9UTHR1NdHQ0ixYtYvPmzR3232dNbiMmJoZNmzYRFRVFcnIyrq6uxMTEEBYWhsFgoLCwkN9+\n+w24/QS7+vp6JEli5MiRmmn56quvWLBggc16du7caVW/YjHEDhBTKDvArk1OSEjobgldQp++JxcV\nFbFjxw7c3Nzw9/fnm2++ISwsjPPnz5OamookSRiNRvbs2YNer2f79u1MnDiR4OBgTfSUlJTw2Wef\nIcsy/v7+lJaWWrRNmzaNpKQkwsPD+f3330lLS+PkyZOcPXuWxsZGpk6dyuDBg8nMzESWZSIiIoiM\njLSq3z59Ju/YsQO9Xo9er+fs2bMMGTKEpUuXEhQURHFxMQAjRoygvr6eq1evcuHCBc0MbtOzfv16\ntmzZwpkzZ9ppM5vN+Pv7s2TJEoKCgigqKmLXrl24uLjg5eVFfn4+NTU13Lx5k4iICMaOHWt1v33a\n5DuRJAk3NzcAHBwc2k2pYmNjSUpK4vnnn9dUw50PAT99+rTl3zqdDlmWLXPhNn06nY6VK1eybNky\nQkJCeOSRR0hMTOTGjRts3brV6n779OU6ISEBg8GAp6cnoaGhlJWVAdyVQz5+/Hjee+89zU2Oj4/H\nYDCg0+lYsWKFRdvTTz/dLo1KkiQkSWLevHkkJSUBt6dftbW17N27Fz8/P0JDQ63u1+6nUC0tLaxa\ntYqpU6fywgsvdLccTbB7k+0Bu7kn2zPCZDtAmGwHCJPtAGGyHfA/e+RONmpdHdgAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb6c36c0d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set(font='Arial',style='white',context='paper',rc={\"xtick.major.size\": 0, \"ytick.major.size\": 4,\"mathtext.default\":'regular'})\n",
    "with matplotlib.rc_context({\"lines.linewidth\": 0.7}):\n",
    "    plt.figure(figsize=(3/2.54,4/2.54))\n",
    "#     plt.figure(figsize=(3,4))\n",
    "    ax=sns.pointplot(x='cone_input',y='area',data=ac_connectivity,ci=95,order=[0,1],\\\n",
    "                     linestyles='',markers='s',scale=1.5,color='black')\n",
    "    ax.set(ylim=(0,27),yticks=[0,5,10,15,20,25],xlabel='',xticklabels=['Rods\\nonly','Rods +\\ncones'],xticks=[-0.1,1.1])\n",
    "    ax.set_ylabel('AII contact area [$\\mu m^2$]')\n",
    "    sns.despine()\n",
    "    ax.spines['left'].set_position(('outward',3))\n",
    "#     plt.savefig('figures/rbc_AII_contacts.svg',bbox_inches='tight',dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Density from flattened OFF- and ON-SAC (Supp. figure 3A, top)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "volume_density=np.loadtxt('data/density_data_BC_flattened.gz')\n",
    "volume_density=volume_density[:,BC_IDs[:,3]==71]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANIAAAByCAYAAAAmhY+sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl0FFX2xz+9Jr2ks3WWzp5OAANhR2EUBRHRA8qIghvg\nhgvoKOooMzqARpFFjnoUBxgBEURBwYOAjvsoiCyyhC1hMSaBrJ2ts3R6T9fvD07Xz8gWoBNIUp9z\n+uSku6vurer61nv13n33ygRBEJCQkLgo5JfaAQmJjoAkJAmJACAJSUIiAEhCkpAIAJKQJCQCgCQk\nCYkAIAlJQiIASEKSkAgAkpAkJAJAwIXk9XopLi7G6/UGetcSp6GsrKzF3y0vL29FTzo3ARdSeXk5\nN9xww2X7oz3xxBMXvO2RI0d49913A+jNxTNz5swWfa+0tJRFixa1sjedF2VbGisrK6OxsfGCt9fp\ndJhMptN+9u6775KXl8egQYOw2+0UFhbS2NjII488AsCiRYuIjo6mqqqq2XarV6+msLCQiooKRowY\nQWRkJEuWLKF///4cPXqUN998k6VLl1JcXIzNZiM1NVXc1uv1MnfuXORyOTabjaysLLKysggODqau\nro6XXnqJ5557jl69elFcXMzw4cPp3r07ixYtQqlUEhQUxLRp08T9bdq0iX379lFZWcmTTz5Jbm4u\nO3bswOv1MmbMGCwWC1u3bsVsNmOxWLjrrrs4fvw4W7ZsoaSkRDyOm266icGDBzNv3jyCgoJQqVSY\nzWYOHTrEb7/9xvvvv09oaCgul4uXXnrpgn8Pif+nzYRUV1fHlClT8Pl8F7wPuVzOihUrCA0NPe3n\nY8eOpX///kybNo0FCxaIF7NKpWL27NnodDomTZrUbJvMzEy6dOlCdnY2P//8M2PGjKFnz55MnjyZ\n6dOnU15ezqFDh3j77bfJzs7ml19+EbfduXMnKSkpTJgwgby8PLZs2UKXLl24//77+f7779mwYQMu\nl4tHH32Uqqoq3nnnHfbu3QucvCkUFhZSWVlJVFQUAF9++SWLFy+mvr6e+vp6Nm7cyLJly/B6vUye\nPJlbbrmFq6++mjvuuINJkybRo0cPkpOTue666zh48KB4HJs3b6a+vp5hw4Zx/fXXk5ubS2hoKDk5\nOSQnJ1NeXk6fPn3IzMy84N9CojltJqTQ0FAWLVp00S3SmUQEYDAYEAQBmUwmvufz+VAoFDQ1NQGg\nVDY/5DfeeIMHH3yQnj17kp+fjyAIaDQaABQKBXa7HblcftptPR6PaKumpgan09nMtn8b/0sQBARB\nYOTIkQwcOJCNGzc2Ox7/c2VTU9MZu8Z+3/w++YP3/3wcHo9H3MZisRAWFoYgCDQ1NTFt2jRqa2uZ\nO3cub7/9NhEREafYcbvd2O127HY7JpPplOOSaE6bdu3O1C0LJFqtlt69e/Pqq6/idDp58sknkcvl\nvPzyyxiNxlOEHBERwY4dO1CpVDQ2Nja7YGQyGRqNhr59+5KVlYXL5SIuLk78fPDgwcycOZNZs2bh\ndruZMWMGWVlZzJs3j8bGRqZNm8ZPP/3UzN6ECRN47bXX+O6771Cr1YwePVr87Oabb+bll1+mrq6O\nxx9/nNGjRzNjxgxkMhmPPvoopaWlpxyvTqdjw4YNpxzH6NGjefXVV9m6dSs6nY4rr7ySY8eOiV27\n2NhYkpOTT7kxNTU1UVBQgM1mQ6FQoFQqqa6uxmw2o9PpLuan6dDIAr0eqbi4mBtuuIEffviBhISE\nQO5aog04duwYp7skHA4HarWajIwMqXU6DW3aIklc3litVlwuF06nkx07dlBTU4PX6+Waa64hNTUV\nn89HSUmJdIM8DZKQJERKSkrYsWMHX3/9NXK5nJiYGOx2OwsXLmTAgAEMHToUr9crCek0SEKSAKCo\nqIiSkhK++OIL+vbty+jRo9HpdPh8PrZv384333zDnj17GDduHLGxseJIo8RJJCFJ4PP5sFgsrF27\nFpPJxJ133olCoQBOjg5ec801XHXVVXzwwQds2rSJ1NRUjEaj9Kz0B6RYuwugo4XllJaWsmXLFqqr\nqxk3bhwejwe1Wk1wcDCCIFBTU0NVVRU333wzbrebrVu3UlhYeKndvqxoUyHV19dTWVl5wa/6+vqL\n9uFiQoT8tHVYzsMPP3zR+zgbv/32G5s3b+aqq64CoFu3bhgMBrxeL2q1moiICEJCQtDr9fTu3Zvs\n7GwOHTokzs1JtGHXzuFw8Mknn5x2aLWlyGQyJkyYIE5K/pELDRG61GE5GzZsOKu9V155BYBPP/0U\nrVbLLbfcwjPPPMPcuXMJCgq64HPpp6qqih9++AGtVsvgwYNJSUmhpKQEhUKBXq8XvxcWFobb7WbQ\noEHs37+f3NxcwsLCuPbaay/ah45AmwlJo9Fw11134XK5LngfQUFBpxWRnwsJEbrUYTnnsgcnbyC3\n3XYbTz/9NN26dSMpKSkgIvL5fBw8eJCDBw8yYsQIwsPDUavVqNXq035fqVRiMpm4+uqr2bZtG717\n92bfvn306dPnon1p77Rp185gMBAVFXXBL4PBcM79n2+I0MWG5TQ2NtKzZ0+AU8JyZDJZs7CcpKQk\n5s6dS01NzVmP43T2/JOh8+bN49577z3r9i3F37L6ozeMRuMZv1taWspXX33F3r17SU1NJSgoiJ9/\n/hmbzUZ+fn5A/GnPdLhRu/MNEbrUYTktsee/MYwaNYrCwkJiYmJafD4aGxuprq7G7XaTlJQktjaV\nlZXs27ePvXv3cv311xMZGUl1dTVqtRq9Xo/T6aSyshKLxYLdbqeuro6YmBi0Wi2FhYUMGzaML7/8\nktTUVDGUKCEhQRR/Z0MKEWon/Pbbb7zxxhtMmzYNs9ncom3cbjfFxcVijFx9fT0hISEUFBRQW1vL\n2rVrUSgUPPDAA+Tn5+NwOICTLbvdbsfr9RIaGkpQUBBGo5H09HRkMhkWi4UdO3Zgt9vZunUrkyZN\nIjExEa1W22LfOhqSkNo5RUVF2Gw2EhISCAkJafZZQUEBWq222XtWq5WqqirWrl3LiRMneOCBB1Cp\nVJSXl3PllVfidrs5duwYwcHB9OvXj+Dg4NPazcvLIycnh+rqaoqLixk/fjzx8fGo1eo2CU6+3Ohw\nXbvORHl5OYIgEBERQVVVFRUVFURERBAeHk5DQwMymYzGxkbq6+vRaDSEhoZSUVHBmjVrqKioYPz4\n8RiNRnJycujevTvR0dEALboBpqen43a7gZPPobt370apVBIbG0t5eTmxsbGteuyXG5KQ2ikOhwOH\nwyG2OP6/1dXVaDQaLBYLDoeDnTt3iospU1JS+Pbbb6murmb8+PFotVry8/MJCQkhLS3tnDbtdjtu\ntxu9Xo9SqSQjI4OmpiZx0raiooLg4GC0Wu0pgz4dnRYJqaCgoNkSa4lLT3l5OVqtlvr6eoqKioiL\niyMsLIyQkBAKCwvRaDT8/PPPREREkJmZyeHDhyksLESr1XLDDTdQWVkJgNFopE+fPmccJHA6nXi9\nXgRBID4+nuDgYEpKSnA6nQQHB9O1a1e8Xq+4jkmn06HRaCgvL+9UXbwWCWnJkiWUlpZyzTXXMGrU\nqGaL2yTaHv/Io9VqZfv27Xg8HvLy8ggODiYuLo7MzEz27NkDQP/+/fF6vYSEhLBt2zZSU1MJDg6m\nS5cuhIaGEhkZCZxsbWQymdh6+Xw+goKCCAsLQ6PRoFQqRbElJCRQWVmJzWZDq9WSlJSEzWYTl+Ur\nFApiYmIkIf2Z2bNn4/V62bZtG7Nnz6aqqopbbrmF22+//ZSHWYnA4Q+LiomJEaMMrFYrNTU1KJVK\ndu7ciU6no0uXLlRUVNDQ0EB+fr44r5OSkkJlZSUNDQ2sW7eOmpoahg0bdsoEqsfjQa/XixHdTqcT\nlUolBq6ejqioKHH1bGRkJKmpqVitVpxOJ0eOHEEmkxEeHi4+d3V0WjToX1xczIoVK1i6dCkajYbJ\nkydjMpl47LHHWtu/TktDQwONjY3o9XosFgtFRUWUlZVRW1uLVqsVk6gkJSXhcDjwer0YDAbi4+OB\nk3Nd1dXVHDx4kJUrV1JWVsbYsWNP+yzkdrubLYsIDg4+q4j8hIeHk5iYSENDA0ajkczMTDEEy2Kx\nkJ+ff1HJbtoTLWqR5s6dy6hRo3jvvfeaDYdardZWc6yzU1VVJbb2/tbI5/Oh1WqxWCxUVlZiNps5\nceIEX375ZbPfQqvVYrfbxf8jIyN5+OGHMRqNhISE0NTUhEwmExO7XExXPSgoSBypMxqNDBs2jK1b\nt4riLCoqIjk5+YL33144q5DcbjdNTU3o9XqGDh2KIAg4HA4xlm3s2LFt5Wenor6+HrlcTl1dHVar\nlZiYGDQaDTabjZqaGo4ePYper2f//v1s3rwZk8nEmDFjkMlk2Gw2ZDIZOp2OtLQ0ZDIZWq2W6Oho\nQkJCqK6uRqVSodFoSEhICEjMnl6vx2w2U1hYSFhYGFdccQVVVVXs2bMHr9dLUlJShx/BO6uQvv/+\ne1avXs2RI0coLS0VhzT9sWUSrYPVasVqtbJ79258Ph/BwcEkJydTW1uLxWIBTj4/7dy5kyFDhnDj\njTcil8sRBAGPx0NVVRUmk0kcho6MjCQ8PBzgrOnMLga5XE5KSgrHjx8nOjoai8WCQqGgqqqKAwcO\n0Lt371axe7lwViGNHDmSkSNHsmfPHvr3799WPnVqSktLyc/Pp6CgAIPBQFxcHHl5eRw9ehSA2tpa\ncnJyqKurY9iwYQwfPpzY2NizBpy2FXK5nNTUVFwuFxkZGRw5coTDhw8TEhKC2Ww+JfKiI3FWIb30\n0ktkZWUxe/bsU5rmdevWtapjnQ2LxUJVVRUFBQWUlZWJk5wHDhwQQ3Z69OjBX/7yF6655hrq6urI\nyMi4LOf3kpOTcTqd6PV6mpqa2Lt3L3q9nmuvvbbDdvHOKqSsrCwAPvvsM/G9zhj+0Zr4fD5yc3Np\naGjgl19+4ddff2XIkCEUFxeTk5OD0WhkwIABZGZmEh0djVarxeVyYTabL9t5Go1GQ0REhLiw0WKx\ncOLECX777Te6du16qd1rFVo0ardw4UJMJhPFxcXs3r2b5ORkceWmxIXj8/k4fPgw//vf//jpp5+Q\ny+Vce+21qNVqRo4cyfDhw5HJZMjlckwmE6GhoeJQ8+U+PxMZGYnD4SA0NJSuXbuyZcsWEhISMBqN\np02R3N5pkZC2bt3Kxx9/zAMPPMCKFSuYOHFia/vV4fFPnm7cuJFdu3YxZMgQ0tPTqaysxGQyERQU\nxBVXXEFERAQOhwOr1YrP5yMtLa1drPnR6XRotVrS09Opr68nIiKCrVu3EhISQnh4eIfr4rVISIIg\nsHz5crp3786BAwdwOp2t7VeHxm63c+DAATZu3MiRI0cYO3YsSUlJ4jJ8s9lMWlqauAhPr9c3y5/Q\nXjCZTLhcLmJjYxEEge+++47o6GiioqJISkq61O4FlBbd2l577TXUajVPPfUUJSUlzJ8/v7X96tBs\n2bKFBQsWcPz4cSZMmED//v0JDw/nxIkTpKeno9Vqz5g3oT2hVCrFJfJyuZx+/fqxZcsWcQSyI9Ei\nIfmXUi9dupTff/+dL774orX96rAcO3aM1atXEx4ezpNPPsnAgQNRKpXk5OSIIT7ns5T8cic2NhaZ\nTEZ0dDTR0dG43W62b9/eLvL9nQ8tEtIzzzxDeHg4mZmZ4kvi/HG73axatYrGxkbuvvtuevTogdvt\n5vfff8fhcNCvXz80Gs0ZV6W2R9RqNaGhoXTv3h25XM7QoUPJzs5m165dl9q1gNKiZ6SwsDDGjRvX\n2r50aARBYN26dezbt48xY8Zw9dVXi3NHxcXF9O3bF6BDLlGJioqivr6eAQMGsGPHDrp3784vv/zC\nlVde2WGmUlokJJVKxSOPPCImvwCa1T6VODuCIPDtt9+yfv16evTowZ133ikuyDt48CCJiYmEh4eT\nlJTULkbkLgS9Xo9GoyEpKQlBEPjxxx/59ddfufXWWzvECF6LhPTggw+2th8dlsbGRn766SdWrlyJ\n0WhkypQpqFQq8vLyyM3NJS4uji5dupCYmNghBhjORHR0NHl5eWRmZlJSUsIVV1zB9u3biY+PJyoq\nCpfLhSAIGAyGdtlKtej2l5SUxKZNm1i1ahWlpaUBiRjuDBQVFbFy5UqWLVtGdHQ0TzzxBJGRkeTn\n53Po0CFxniUqKqrDn1OZTCYu4UhNTcVkMlFQUCBeTwaDgdDQUJxOJ3V1dZfa3fOmRUL617/+xb33\n3kttbS39+/dn1qxZre1Xuyc/P5+srCy++uorevTowbPPPktCQgLFxcUcPXoUl8tFZmYmqamprRaR\nfbkRGxuLXC4nPT0dlUpFZmYm33//fbN5SY1Gc85MtJcjLRKS0+kkIyMDgMTExA41qtQa2O125s+f\nj91u57HHHmPixIk0NTVRWlrK7t27qaqqIjMzk/j4+A7fEv2ZhIQEPB4P6enpREdHU11dTW5uLj6f\nD4fDgSAIKBQKqqurL7Wr50WLnpEGDhzI888/j8ViISsri169erW2X+0WQRD4z3/+Q1lZGVOnTiUj\nI4PCwkJKS0vFVFmDBw8mJCSkQ8acnQt/LgeFQkFhYSGDBw8mPz+f4uJiBEEgPDycQYMGYbVaUalU\n58z3frlwzhbpk08+ITc3l82bN3P8+HFiYmJ4/vnn28K3docgCCxatIgff/yRYcOGkZKSwsGDBzlw\n4ACCINCzZ0+uuuoqwsLCOuQwd0uJiIjA6/Vy/fXX09TURElJCTExMfTr1w+bzUZ2djZ6vZ7q6mqK\niooutbst4qxCmjt3LpWVlbzyyits375dLF0ihQidSlNTE7Nnz+brr79m8ODBXHfddRw4cICCggK6\ndu1K7969iYuLIzk5uUNFLlwo/iSTV111FceOHePQoUPEx8fTt29fysvLyc/PR6vVIpfLKSgoEKuG\nXK6ctWu3b98+1qxZI/4fExPDtGnTpOjvP+FwOMjKyuLw4cOMGDGC7t27s2vXLtRqNb169cJsNl/2\nyx7ampiYGH7//Xfi4+Pp06cP2dnZDB8+nOTkZFJTU8nJyUGlUpGYmIhGo6GgoACVSnXZzrWdVUhn\nmtfoCBNogUAQBL7++mtWr16NzWbjtttuIygoiLy8PMxmM/Hx8aSkpJy1OFpnRSaTERUVRWVlJVdf\nfTV79+7lq6++4oEHHiAzMxO32012drYY9eGvfVVYWIhSqRSvQZ1Od1k8a55VSBaLhY8++ui073d2\nKioqeOutt7BarXTp0oWUlBQcDgd6vZ5+/foRHx8vBmxKnJ7Q0FBCQkLQ6XQMGjRIrALYr18/+vfv\nT3x8PPv37+fHH3/EbDaTlJQklqjx09DQQENDA4mJiZe0pTqrkCZPnnxe73cGioqKWLVqFQcPHmTg\nwIFi9TqNRkNaWhpJSUnExMSgUqkutavtAv/q37/+9a/k5ubyww8/kJ6eTnh4OCaTicjISHJzc8UE\nMGlpaSQmJqLX61EoFGLl9fz8/GbZktqaTl8fyZ+rTy6Xo1Kp8Hq9VFdXk52dTV1dHXK5HLvdzvHj\nx8Wkjf65tKCgIPr374/BYEClUkkt0EXg8/n4/PPP+eCDD7j99tsZOnQoMplMnLP0eDzk5+dz9OhR\nBEFApVIRFRVFQkKCeN5dLhcej0fs+ikUCiIiItqka93uy7q4XC5yc3MpLCzkxIkTuN1ufD4fMTEx\nYn616upqnE4nCoUChUKB2+3GbrcjCAI6nY6wsDAEQcBms+F0OmlqasLj8YiTpUFBQZhMJjEdcFxc\nHJGRkaSlpUmDCAFCLpdz4403kp2dzX//+18GDBhAYmKimEFWo9GIhagbGxupqKigsrKSX3/9ldjY\nWHr16oVGozllgru8vJympiaUSmWrDlRcdi1SQ0MDJ06cwOfzNStm7L8L2Ww2CgsLKS8vp7a2FqfT\nKYaY+DOI+gsuq9VqwsLCTnvy/piTWqFQNKvE4BecRqMR7RqNRsLDw0lNTe3Q+dkuNTk5OcyePRuj\n0cgLL7xAbGwsDQ0N1NbWijc3/+8rCAJlZWUcOHCApqYmevTocUaxCIJAQ0MDZrO5RXnNz5dWE9Kq\nVatQq9U4HA7q6urEKt+xsbE4nU4cDgdKpRKlUonP5xMry/mz5vwRj8fT7MJXqVSnPVkymQyVSiUu\ncVYqlURGRoqL5VQqFSEhIRiNRinM6TJmw4YNLF++nG7dujFlyhRSUlLEz/wF1OD/q7+73W4OHTpE\nUVERCoUCvV6PVqslODiY8PBw8RqQyWQ0NDSQkpJySnX7i6XVhHTHHXdgMBjEVLr+MiFNTU3NROC/\n4OFkjFp9fT0ej0cMoQkPDycoKEhsMZRKJTqdDp1O12wNj8fjEVsS6TmlfeNyuVi/fj2ffvop0dHR\nPP7446eEpdXV1VFTU4NcLhdvinV1dVRUVGCz2bDb7TgcDrGWVFBQEElJSaSlpeFyuUhOTg7ospVW\nE9LUqVOJjY0VW4+goCDUajVOp1MUjkKhQKvVEhUVRWxsrNRKSIi43W5+/PFHVq1ahcfj4ZFHHuGG\nG2445Xv+GlJnyrLk7xFVVFRw4sQJ4GSWptjY2IDO8V12z0gSEn/k6NGjLFy4kOPHjzNkyBBuuukm\n0tPTm7UmLpeLsrIyfD7fKfNMf8TlcpGXl0d+fj4qlYr09HQSExMxGAyo1Wqx+3chSEKSuOyprq5m\n2bJl7NixA6/Xi16vJzo6mri4OEwmE2azWUzS39DQgMfjQRAEsTqHv0ekUqnE6Y5Dhw5RVlaGVqvF\nZDKJjx5BQUEMGDDgvH2UhCTRLvD5fBw9epQDBw5QWlpKWVkZFRUVWK1W/JewXC4nLCwMo9Eo1rBN\nTk4mMjISvV6PWq0Wn9M9Hg8Wi0VMQONPzgnw6KOPnrd/7X4eSaJzIJfLycjIEBeYAuLcX0FBAeXl\n5WIlw8rKSo4dO8aOHTvweDzN9qNQKJDL5cjlcnFgyv+/Wq2+4AEISUgS7RZ/HohevXqddrGpf+7I\nL7Da2lqampqavbxeLz6fD6/XK753IUhCkuiwyGQyDAYDBoOBLl26tKqty29hh4REO0QSkoREAJCE\nJCERACQhSUgEgIAPNvhHPTpa2Q6JzkVsbOx5BbYGXEiVlZUAjB8/PtC7lpBoM843oCDgQsrMzCQl\nJYX33nuvVdZ9nInJkyezePHiNrN3Ke12FpuXyu7kyZPPO5F/wIUUHByMTqcjOTk50Ls+K2q1+pKE\nJF0Ku53F5qWy+8elPS1FGmyQkAgAipdffvnl1tjxpSiPealKcnaWY5XO75kJePS3hERnROraSUgE\nAElIEhIBQBKShEQAkIQkIREAAjKPVFpayuOPP05GRgbR0dEYDAZKSkqw2Wy88MILrZKPuaSkhIUL\nF6LX6zEYDAQHB7e6TYCPPvqIQ4cO4fF42Lt3L+PHj291u0eOHGHx4sVial6j0djqNvfv388HH3xA\nZGQkvXv3pqKiotVtHj9+nKeffpr169ezbNmyZvbcbjfz5s0jLCyMtLS0gEXO/NGmz+dj6dKlyOVy\nHn74YSwWS4ttBqRF2rVrF1FRUQD06dOH3bt3M3PmTO644w4+/fTTQJg4heXLl5OUlER9fT09e/Zs\nE5twMvRpzpw5mEwm3nzzzTaxGxERIS6lDg0NbRObX3zxBVOmTGH69Ol88803bNu2rVVtVlVVsW7d\nOrRaLW63m127djFz5kzGjh3Lp59+yieffMJ9993HzJkz2bx58wWvZD2TTThZndLpdIo5IM7HZkCE\n1KtXL+bMmcOcOXPEuxicDPzzx94FGn96plmzZrF06dI2seknPz8fr9cr5gBvbbtr1qxh6tSpvPXW\nW/z6669ia9CaNu+//36WL18uFpX2VxlsLZtGo5G///3vaLVa6urqMBqNwMmCZBUVFVRVVWEymQAw\nGAw0NDQE1CbAPffcw6BBg8TPz8dmQIR0+PBh3G43cDJEqKKiAjgZAd5aSeajoqLQ6XRiaY/a2tpW\nt+nn448/ZuLEiURGRraJXZfLJYpHr9dTWlra6jbLy8t54okneP755/H5fGKV8bY4v38+rzExMcTF\nxVFWVgZAfX19QIs0n2kqNS4uTlzFcC6bAXlGSk5O5vXXX8doNDJ8+HCcTidZWVnU19eTlZUVCBOn\nMGnSJN544w1CQkIYPXo0Vqu11W36+f3338ViygMHDmx1uxMmTGD+/PlERETQp08fVCpVq9uMiYlh\n3rx56PV6Ro0a1Sa/qR+5XH7KeXU6ncyZM4f169dz4403BrSqxJmSQo4dO5Y5c+ag1+vPaVOKbJCQ\nCADS8LeERACQhCQhEQAkIUlIBABJSBISAUASkoREAJCEdJGUlJTw1FNPATBixAjuu+8+7rvvPu6+\n+25++OGHZp+fi9tuu+28bH/++ecAvPvuu2zevLnF2/3zn/9k3LhxVFdX89hjj3H//fdjs9kAmDFj\nhjgnOH36dIYMGXJePnVWpNzfAcA/D2EwGFi5ciUAjY2N3HfffSxYsKDFxavOt8jVhx9+eN7i89uZ\nM2cOZWVlDB48mPj4eLZt20Z8fDxdunQRKzLMmjWLsWPHnvf+OyOSkALIH6fk6urq0Gq1Z5w197N2\n7VrWrFlDamoqXq8XOJkKasmSJcjlcqZOncrAgQO59dZbSUxMpLy8nGeffZbGxkYKCgp46623UKlU\nrF69mmXLlhESEsI777zDN998w8qVK/H5fPztb3/juuuuO8W2RqPB5XLhcDgwGAysWbOGGTNmBPak\ndBKkrl0AsdlsYtfulVde4cUXXzznNh9//DFr167l2Wefxe12IwgCixYt4sMPP+T9999nwYIFAFRU\nVPD666+zfPlyFi9ezE033YTZbOaZZ54BoF+/fqxcuRKDwUBOTg5ff/0106dPZ+nSpaJA/0xaWhq1\ntbVkZ2cTGhpKWloaq1ev5rXXXhO7ehItQ2qRAkhISIjYtfNTUlJyxu9XVVURGxuLXC4nLi6O8PBw\nampqKC4uZtKkSQiCgNVqxe12k5qaKhYc9leX+2Nr161bNwDCwsJwuVw899xz/Pvf/6ayspI777zz\njD4899z1cG4iAAABZ0lEQVRzwMlno6lTp7JkyRLGjBnDpk2buOeee87ZokqcRBJSADnTRXem98PC\nwigvL8fj8VBTU4PVaiUiIgKz2cyKFSvwer0sX74ctVrN8ePHcTqd2O12MfHmmZ6pBEFg48aNzJgx\nA7VazdixY7n55pvP6Hdubi5ms1nM5SaTyZqVgpQ4N5KQAsiZLuydO3eKD+0ZGRm8+uqrACiVSh56\n6CHuuusuUlJSCA0NRSaTMXHiRMaPH4/T6eTuu+8Wvzt16lRqamr4xz/+AZwMFn7ttdcIDQ1t5oNM\nJqNbt25MmDABvV7PxIkTz+r3mjVrePHFF8Uo+lmzZjF//vyzHpNEc6Sg1XbCrbfeyqZNmwKyrxde\neIGHHnqoRVXs7rjjDj777LOA2O3ISIMNnZQXX3wRq9V61u9Mnz5dXIckcXakFklCIgBILZKERACQ\nhCQhEQAkIUlIBABJSBISAUASkoREAJCEJCERAP4P8LxI4Ts482sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb6bb513b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8/2.54,3/2.54))\n",
    "sns.set(font='Arial',style='white',context='paper',rc={\"xtick.major.size\": 4, \"ytick.major.size\": 0})\n",
    "plt.fill_between(x_scale,(np.mean(volume_density[:,rbc_cone_input-697],axis=1)-np.std(volume_density[:,rbc_cone_input-697],axis=1)/np.sqrt(len(rbc_cone_input)-1))\\\n",
    "                 /np.sum(np.mean(volume_density[:,rbc_cone_input-697],axis=1)),\\\n",
    "                 (np.mean(volume_density[:,rbc_cone_input-697],axis=1)+np.std(volume_density[:,rbc_cone_input-697],axis=1)/np.sqrt(len(rbc_cone_input)-1))\\\n",
    "                 /np.sum(np.mean(volume_density[:,rbc_cone_input-697],axis=1)),facecolor='0.5',alpha=0.3)\n",
    "plt.fill_between(x_scale,(np.mean(volume_density[:,rbc_rod_input-697],axis=1)-np.std(volume_density[:,rbc_rod_input-697],axis=1)/np.sqrt(len(rbc_rod_input)-1))\\\n",
    "                 /np.sum(np.mean(volume_density[:,rbc_rod_input-697],axis=1)),\\\n",
    "                 (np.mean(volume_density[:,rbc_rod_input-697],axis=1)+np.std(volume_density[:,rbc_rod_input-697],axis=1)/np.sqrt(len(rbc_rod_input)-1))\\\n",
    "                 /np.sum(np.mean(volume_density[:,rbc_rod_input-697],axis=1)),facecolor='0.7',alpha=0.3)\n",
    "\n",
    "plt.plot(x_scale,np.sum(volume_density[:,rbc_cone_input-697],axis=1)/np.sum(volume_density[:,rbc_cone_input-697]),label='rod and cone contacts',c='0.3')\n",
    "plt.plot(x_scale,np.sum(volume_density[:,rbc_rod_input-697],axis=1)/np.sum(volume_density[:,rbc_rod_input-697]),label='rod contacts only',c='0.6')\n",
    "plt.xlabel('IPL depths [%]')\n",
    "plt.ylabel('Density')\n",
    "plt.yticks([])\n",
    "plt.xlim(50,110)\n",
    "sns.despine()\n",
    "plt.legend(bbox_to_anchor=(0.65, 1.13))\n",
    "# plt.savefig('figures/rbc_density_comparision.svg',bbox_inches='tight',dpi=300,transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Density aligned on ON-SAC (Supp. figure 3A, bottom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rbc_densities=np.loadtxt('data/density_data_RBC_ON_SAC_alignment.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANIAAAByCAYAAAAmhY+sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl0U9X6v5+MTdJMnZu0tKUDtUgBqYjiiKB4UbkgIA7g\nPOCIepWrfgFFQUCXuhyucBVF9Io4XRW4Lr2KgqBlkslSqJ0oNG06pUOSNnN+f3T1/OyFQoG0pe15\n1so/yTnnfc9JPtl7v/vd75YEg8EgIiIip4W0px0QEekLiEISEQkBopBEREKAKCQRkRAgCklEJASI\nQhIRCQGikEREQoAoJBGRECAKSUQkBIRcSD6fj/Lycnw+X6gvLXIMKisrO32s1WrtQk/6NyEXktVq\nZezYsWfsl/bAAw+c8rkHDx7kzTffDKE3p8/8+fM7dVxFRQXLli3rYm/6L/LuNFZZWYnT6Tzl88PD\nwzGZTMf87M0336SoqIjzzz+f5uZmDh06hNPp5O677wZg2bJlxMbGUltb2+68jz/+mEOHDlFdXc2V\nV15JVFQU77zzDjk5ORQUFPDKK6+wYsUKysvLcTgcDBw4UDjX5/OxZMkSpFIpDoeDBQsWsGDBAlQq\nFY2NjTzzzDM8/vjjDB06lPLycsaNG8fgwYNZtmwZcrmcsLAw5syZI1xv3bp17Nmzh5qaGh566CHy\n8/PZunUrPp+PyZMnU1VVxZYtW0hNTaWqqorp06dTVlbGzz//jMViEe5j/PjxXHTRRSxdupSwsDAU\nCgWpqank5eVRWFjIe++9h8FgwO1288wzz5zy9yHy/+k2ITU2NnLfffcRCARO+RpSqZRVq1ZhMBiO\n+fnUqVPJyclhzpw5vPHGG8KPWaFQ8MILLxAeHs6dd97Z7pwhQ4aQkZHB7t272bx5M5MnTyY7O5tZ\ns2Yxd+5crFYreXl5vPbaa+zevZtffvlFOHfbtm2kpKQwY8YMioqK+Pnnn8nIyODWW2/lhx9+4Ouv\nv8btdnPPPfdQW1vL66+/zq5du4DWP4VDhw5RU1NDTEwMAP/5z39Yvnw5TU1NNDU1sXbtWt599118\nPh+zZs3immuuYfTo0UyZMoU777yTs88+m+TkZC655BJ+//134T42bdpEU1MTl19+OWPGjCE/Px+D\nwcD+/ftJTk7GarUyfPhwhgwZcsrfhUh7uk1IBoOBZcuWnXaL1JGIAPR6PcFgEIlEIrwXCASQyWT4\n/X4A5PL2t/zyyy9z++23k52dTUlJCcFgELVaDYBMJqO5uRmpVHrMc71er2DLZrPhcrna2W47p+0V\nDAYJBoNMmDCBUaNGsXbt2nb30zau9Pv9HXaN23xr86ktef9/78Pr9QrnVFVVYTQaCQaD+P1+5syZ\nQ0NDA0uWLOG1114jMjKyw2faEYFAgLq6OoxGIwqF4qTP72t0a9euo25ZKNFoNAwbNoznn38el8vF\nQw89hFQq5dlnnyU6OvooIUdGRrJ161YUCgVOp7OdECQSCWq1mnPOOYcFCxbgdrsxm83C5xdddBHz\n589n4cKFeDwe5s2bx4IFC1i6dClOp5M5c+awcePGdvZmzJjBokWL+P7771EqlUycOFH47KqrruLZ\nZ5+lsbGR+++/n4kTJzJv3jwkEgn33HMPFRUVR91veHg4X3/99VH3MXHiRJ5//nm2bNlCeHg4I0eO\n5I8//hC6dvHx8SQnJx/3j6kjqqqqqKysRKFQYLFYiI6OJjEx8aSv05eQhHo9Unl5OWPHjmXDhg39\n/uH2FWpqarDZbKhUKux2O0C7Vsjn86HVaklOTu4pF3ucbm2RRHofDQ0NOBwOjEYjFouFhoYGvF4v\n+/fvx2azkZmZyYgRI6ivrychIeGo7m9/oX/etUin8Hg81NXVodVqBQEVFxezadMmQVz79u1j8+bN\nzJw5k8rKSgYMGNDTbvcIopBEjonb7ebw4cPo9Xrcbjc2m421a9eyf/9+MjMzufLKK9Hr9VgsFtav\nX8/y5ct57LHH+q2QxBQhkaMIBoMcOXIEvV5Pc3Mz5eXl7Nixg+LiYqZMmcINN9xAZGQkcrmcqKgo\nbrvtNlpaWti4cSMul6un3e8RRCGdAn09LcdqtaLRaLDb7ZSXl5OXl4dKpWLcuHFoNBrKysqwWq1I\nJBI0Gg1er5dzzz2XnTt3cvjw4Z52v0foViE1NTVRU1Nzyq+mpqbT9uF0UoTa6O60nLvuuuu0r9FZ\n6uvrcblcNDU1UVlZSX5+Pi6Xi9raWmH6wGazUVlZyYEDBzCZTMTHxzNgwACcTid79uxpN4fVX+i2\nMVJLSwuffPIJpxNtl0gkzJgxQ5iU/DOnmiLU02k5X3/99XHtPffccwB8+umnaDQarrnmGh599FGW\nLFlCWFjYKT/LY1FbW4vdbqepqYmGhgYOHDiAx+OhpaWFsWPHotPpiIqKAsBut7Nx40YKCwsZPHgw\n2dnZ/PDDD+zdu5fzzjuPlJSUkPp2ptNtQlKr1UyfPh23233K1wgLCzumiNo4lRShnk7LOZE9aP0D\nmTRpEo888giZmZkkJSWFXESBQID6+nrsdjuNjY0UFBTgdDqRyWRcfvnlwnioDZ1OR3p6OkVFRSQl\nJREbG0t2djZbtmzhyJEj/U5I3dq10+v1xMTEnPJLr9ef8PonmyJ0umk5TqeT7OxsgKPSciQSSbu0\nnKSkJJYsWYLNZjvufRzLnlKpJCsri6VLl3LTTTcd9/yTxePxUFpaSnNzM3a7nd9//x2Hw4FOp2P0\n6NHI5fJjZqVkZGSgVCr57bff8Pv9jB49GplMxtatW4WJ2/5Cnwt/n2yKUE+n5XTGXtsfw9VXX82h\nQ4eIi4sL2fPyeDyUlZWhUqkoKyujoKAAnU5HXFwcycnJ6PV6wsLCcDgc+Hw+9Ho9Ho8HmUyGUqlk\n5MiRbNu2jS1btjB06FCys7PZuXMnFRUVZGZmhszPMx0xRaiXUFhYyMsvv8z999+PyWTCYDCgVCpP\n65rBYJCioiL0ej1lZWVs2bIFtVqN2WwmPj6emJgYgsEgsbGx6HQ6gsEgtbW1KJVKvF4vdrsdjUZD\nU1MTmzdvJioqCrlczj//+U8mTpzI7bffjkwmC9ETOLMRw9+9hIyMDJ588kkMBgM+n4+ysrLTGm96\nvV6KiorQ6XS4XC62b99OeHg4sbGxJCYmYjabCQQCgoigtWWMiYnBYDAQHR0tZLTr9XpycnKoqqpC\nKpWSlpZGbm4uZWVlobr9Mx5RSL2EtmxrlUqFVCpFr9e36waeTMciEAhw6NAhYUy5fft2FAoFzc3N\nDBw4EJPJhN/vJy0tTRDRsTCbzTgcDgDi4+NJT0/HYrFw7rnnUlNTw+7du/F4PKd+070IUUi9ALvd\njsvlOipQolAoKC4uprS0lIKCAmpqajp1vcOHDwtdtR07dmCz2SgtLWXEiBGYTCZaWlpISUk5am3V\n/yKTyTAajTQ3NwMwaNAg5HI5KpWKmJgYfvnll2OO+foinRJSaWlpV/sh0gEej4eqqiohg6C4uJh9\n+/ZRXFyMTCZDq9Wi0WiIiIjA4XCcMOuioqICmUxGMBhk7969WK1W9u3bR3p6OsnJybS0tJCamtpp\n/6Kjo4mJicHpdKJQKMjMzKSuro6RI0dSVFREXl6eEDHty3QqavfOO+9QUVHBhRdeyNVXX91ucZtI\n12KxWNDpdDQ3N/Prr7/idDoJCwvD7XZjtVq54IILkEqlBINBNBqNMH7637VBDoeD2tpa/H4/TqeT\n/fv3Y7fbKS0tRS6XM3LkSILB4EmJqA2tVksgEKC6upqUlBRKS0sJBAJER0fz7bffMnjw4FO6bm+i\n01E7n8/Hr7/+yqeffkptbS3XXHMN1113HRqNpt1xYtQudFitVrxeL1KplM2bN+P1ehk9ejRarZaa\nmhpyc3PR6XSoVCpsNhtJSUkMGTIEj8eDQqEQ5n4cDgdVVVWEhYVRVFREWVkZPp+PwsJC/vjjD26/\n/XYyMjJISUlBpVKdsr9utxuLxYLT6SQ3NxepVMratWuZNm0akyZNOu54q7fTqa5deXk5q1atYsWK\nFajVambNmoXJZOLee+/tav/6LR6PB4fDgUKhYNeuXTidTkaNGoVWqwUgJiaG888/H41GQzAYJCws\njJKSEn777TekUilOp1PoUlVVVREMBvn9998pKCigqamJn376idLSUqZNm0ZKSgoZGRmnJSJozTwJ\nCwsjMjKSqKgolEolUVFR/Pbbb+Tl5Z1WetiZTqe6dkuWLOHqq6/m7bffbvew6+vru8yx/k5NTQ06\nnY6DBw9SWVnJyJEjj6qvEBsbi1arFSJndrudAwcO0NzczPDhw8nNzcXtdiOTySgpKcHn81FXV0d9\nfT3Dhw8nJycHuVzOkCFDThhY6Cxms5ni4mIGDx7M5s2bGT16NOvWraOoqIiYmBjS09NDYudM47hC\n8ng8+P1+tFotl112GcFgkJaWFiGXberUqd3lZ7/D4/FQXV1NQUEBmZmZmM1m/H4/zc3NwvxNdXU1\ngUBAWAOk0WhISkqitLSU7du3I5VKKS0txWAwCMvDx40bR1JSEna7Hb1ez4gRI0ImIkBYWgGtc1+F\nhYWMGjWK7777DoPBQHx8vNCq9iWOK6QffviBjz/+mIMHD1JRUSHksbXllol0DS0tLTidTnbt2oXJ\nZCIzMxOfz4fX6yU9PZ26ujoqKiooKChg06ZNQthbpVKh1WpRq9UMGzaMP/74A51Oh0wmo7a2luuv\nv16I8J1zzjmn3ZXrCLPZTElJCVlZWUIdw8rKSvbu3UtERATnnXdeyJNue5rjCmnChAlMmDCB3377\njZycnO7yqd9TW1tLSUkJSqWSc845B2gtsDl48GBKS0spLy9n/fr1HDhwgPT0dMaPH49cLqekpITG\nxkZiYmKQSqUMGjQIaP1hp6WlYTKZSEtLO24GfaiIjIyksbGRIUOGUFtbS3Z2Nrm5uZxzzjmCyPoS\nxxXSM888w4IFC3jhhReOav4///zzLnWsv3LkyBEOHjyI1WolOTmZhoYGLBYLWq2WXbt2sXPnTr79\n9ls0Gg0zZ85k0KBBqNVqJBIJmZmZeDwe6uvr0ev11NXVkZCQQFhYWJcsvTgeRqMRm81GMBhk5MiR\nbN68mWHDhrFhwwaio6OxWq3Ex8d3mz9dzUknrZ7oAYjh71MjGAxy6NAhduzYgcViEULaarWaw4cP\nU1JSgkajobm5mdGjR3P55ZcTERFBamrqUYmhTqeTyspKdDodXq8XhUIR0ozxzhIIBCgpKUGlUlFd\nXc3OnTtpaGggJiaGCy64ALPZfMKlMb2FTkXt3nrrLUwmE+Xl5ezcuZPk5GRh5abI6dMWmv7xxx/Z\nunUrl1xyCTabjQMHDuBwOIiPj+cvf/kLLpeLxMREzj77bBITEzv8EYaHh5OYmEhJSQlarbZHRASt\n66nS09MpKyvDbDaTkZHBgQMHqKqqor6+HpVK1b+EtGXLFlavXs1tt93GqlWrmDlzZlf71a8oKCjg\n/fffp6SkhCuvvBK1Wk1WVhbjxo1DoVCgVCoJDw/H5XJhNBo7VfpZpVKRlZUV0ojcqdIm6rS0NOrq\n6pDL5Wzfvp0xY8bg9/v7xFKLTgkpGAyycuVKBg8ezL59+/ptyaWuwGKx8OGHH1JWVsY111xDMBgU\nJlqTk5MJBALCStmT5UwQEbQmtw4YMICKigqysrLYsGEDUqkUm81GVVVVn0g569Q3tGjRIpRKJQ8/\n/DAWi4WXXnqpq/3qFzQ0NPDBBx+Ql5fHmDFjgNYMapVKRXR0NMApi+hMQ6VSER4ejkajYdCgQchk\nMvbu3dtnlqR36ltqW0q9YsUKiouLWb9+fVf71S/45ptvyM3N5cILLyQ1NRWJREJcXBxarZbw8PCe\ndi/kxMXFodfrGTBgAG63G5VKhcViOWENi95Ap4T06KOPEhERwZAhQ4SXyOmxd+9evvrqKzIyMsjJ\nyaG8vFxYVNe28VhfJCoqivDwcIYOHYrP5+PgwYPYbLbT2oDuTKBTYySj0ci0adO62pd+QTAYZN++\nfSxbtoywsDAuv/xyiouLSUlJITU19ZT2K+ptxMfH09zczO7du1Eqlfj9fioqKnr1dEmnhKRQKLj7\n7rtJT08XBrB/3vtUpHMEAgE2btzIqlWr8Hq9XHfddVitVhISEhg4cGCfGHR3BpVKRUREBImJiTQ1\nNbF9+3YuvvhiGhsbe+0fSaeEdPvtt3e1H32eQCDA+vXr+eijj9DpdNx22200NjYil8sZOHAgCQkJ\np10VqDdhMpk466yzhMKacrmc2traXiukTo2RkpKSWLduHf/617+oqKjocwmHXU0gEOCLL75g5cqV\nmEwmZs+ejcFgEFaUZmRkHLVAsj8QHR2NRqNBp9Oxd+9e1Go1R44c6Wm3TolOCen//u//uOmmm2ho\naCAnJ4eFCxd2tV99iq+//prVq1eTmprKTTfdhM/nIz8/H51OR3Z2dr/dzDgyMpIhQ4Zgt9spKysj\nEAgQCAR65Q4enRKSy+USsnUHDBjQZen3fZH//ve/fPjhhwwcOJC//vWvFBQUsGPHDgCGDRt2SjuK\n9xU0Gg3R0dEEAgEkEgnbtm1DpVIJVV17E50S0qhRo3jiiSeoqqpiwYIFDB06tKv96hNs27aNt99+\nm6SkJCZNmkRRURFKpZKhQ4dy3nnncdZZZ/W0iz2ORqNh1KhR7N+/n7q6Oo4cOYJOp+t1rdIJgw2f\nfPIJ+fn57Nq1i6amJiZPnsysWbO6w7deTWFhIa+++iqRkZFcd911FBUVodVqGTFiBJGRkf26Jfoz\nsbGxNDY2Eh4eTnV1NXv27CEyMpJAINCr8vCO2yItWbKEmpoannvuOXJzc4UIi5gidHy2b9/OvHnz\nMJvNTJgwQShMP2LECNLT00UR/QmFQoHBYGDMmDHs3bsXr9dLbm4uarUaq9VKMBikpKTkjJ+wPa6Q\n9uzZw4MPPkhsbCwymYy4uDjmzJnDvn37usu/XoXb7WblypW8/PLLjBo1irPOOgun00lycjIjR47s\ns4U/Thez2YzJZCInJ4etW7cKO/95PB6Ki4uFNVlnMsft2nU0r3GmZBWfCbQtyCstLWXz5s3Y7XbG\njBmDSqVi2LBhGI1GDAaD2AodB4lEQlZWFjU1NRQWFlJcXAxAeno6RqMRaN3X6tChQ0RFRZ2R9fGO\nK6Sqqio++uijY74v0sqPP/5IYWEhUqkUk8mEyWQSMhXS09P71STr6aBUKsnOzqa5uZk1a9aQmJjI\nnj17uOyyy4DWLqBCoaCuro7a2lpUKhUKhULY3C0iIqJHo8nHFVJHQYX+Hmzw+/1s2LCB3NxcTCYT\n1dXV6HQ6xowZg9FoRKvVCssgRDpPXFwcWVlZTJw4kZ9++omcnBx+/PFHXC4XmZmZ7Qq3BAIBYcMz\naK1prtFoeqwOxHGFNHny5O7yo8dwu924XC70en27LqvD4aC6uhqJRILb7ebw4cPU1dVRV1dHY2Oj\n8KX5fD5uu+02kpKSevAu+g4ZGRm4XC4sFgsFBQWcffbZxMbGkpeXh9frJTMzE4lEctQ6rbYVxNXV\n1cTGxna7331m60ubzUZhYSEtLS1CsQ+DwYDdbqe8vJza2lqam5vx+XzCTt2BQAClUolCocDlcmG3\n23E4HMI15HI5Go2GyMhIIfsgEAig0+mIjIwkISGB8847r0+uHepJsrOzsdls/Pvf/+b777/n7rvv\nJjMzk4KCAsrLy8nMzCQhIeEoMalUKlpaWqiqqur2OhVnpJDsdjtFRUVIpVJkMhmBQICYmBgkEglO\np5OamhoqKys5dOgQNpsNp9OJy+XC7/ejUqlQq9VoNBrCwsLQ6/Xt5iICgQAOh0Modmm32/F4PGg0\nGoxGI9HR0ULLJJfLkclkhIWFERUVRUJCAsnJycjlcjHfsIvJycmhsbGR1atXs3r1aqZPn865557L\nkSNH2LVrF3l5eUREROD3+5HL5RiNRgYOHIhKpaKpqYnY2NhuDYp12R6y33//PWFhYVRUVFBXV4dG\no0GtVqNSqZDL5fj9fnw+H263G6fTidPpxGazYbfbO1xe3dLSQjAYFHatO+YNSSQoFApBgGq1GqPR\niF6vJzIyErPZjEaj6fB8l8uFx+NBp9OJ0ckepry8nC1btrBu3ToCgQDjx4/n4osvxu/3c/jwYRwO\nB3K5HK/Xi81mQyqVChWWlEplt45Tu0xIU6ZMEX6Mcrlc+Hdv69+2vf6Mx+Ohrq4Ot9uNXC7HYDCQ\nkJCAXq+npaUFn8+HQqFAq9USGRlJYmIiERERKBQKPB4PHo+H8PDwXjMbLnJiHA4H27ZtY+3atZSW\nljJ8+HCmT59+lEhaWlrIz8+nvLyc9PR00tLSjtojqivpsq5dTEwMcXFxhIWFER4ejlKpFLJ7g8Gg\nIC61Wo1Op8NgMGA0Gk95OYFKpRKTafsgbRs4DBgwgG+++YZNmzZRUVHBpEmTGD58uDB2VavVjBgx\nAr1eT35+Pj6fr28Iadq0ab166bDImYNMJmPQoEEMGDCAzMxMPvnkE95++22Sk5O56KKLuOCCCwgL\nC0MikZCRkYHP56OoqIiGhgZhQrer6bKunViyWKSrKC8v54cffuCXX36hqqoKg8EgZNS3VWP6/vvv\nSUxM5IorrugWn0QhifRaGhsb2blzJxs3biQ/Px+v10tsbCzZ2dkkJiZSX19PWloaRqORmJgY1Go1\nMplMqF4rl8tDFlA6I8PfIiKdwWAwMHbsWMaOHUtjYyM///wzO3fuZPPmzQSDQYYPH05LS0uHUxXB\nYBCfz4ff72+XXf7II4+ctC+ikET6BAaDgWuvvZZrr70Wr9dLXl4edXV1wtRK22T8/yKTyZDLT18G\nopBE+hwKhULYoK276BuFpUVEehhRSCIiIUAUkohICBCFJCISAkIebPD7/QC9rpySiMifiY+PP6lo\nXsiFVFNTA8DNN98c6kuLiHQbJ5tQEHIhDRkyhJSUFN5+++1uzcKeNWsWy5cv7zZ7PWm3v9jsKbuz\nZs066SXrIRdS2xaH3Zl5C63FM3oiJakn7PYXmz1lty196GQQgw0iIiFA9uyzzz7bFRfuie0xe2pL\nzv5yr+Lz7ZiQZ3+LiPRHxK6diEgIEIUkIhICRCGJiIQAUUgiIiEgJPNIFRUV3H///WRlZREbG4te\nr8diseBwOHjqqaeIiIgIhZl2WCwW3nrrLbRaLXq9HpVK1eU2AT766COhfO6uXbu4+eabu9zuwYMH\nWb58OfHx8UgkEqKjo7vc5t69e3n//feJiopi2LBhVFdXd7nNsrIyHnnkEb788kvefffddvY8Hg9L\nly7FaDSSlpYWssyZP9sMBAKsWLECqVTKXXfdRVVVVadthqRF2rFjBzExMQAMHz6cnTt3Mn/+fKZM\nmcKnn34aChNHsXLlSpKSkmhqaiI7O7tbbEJr6tPixYsxmUy88sor3WI3MjISq9UqFProDpvr16/n\nvvvuY+7cuXz33Xf8+uuvXWqztraWzz//HI1Gg8fjYceOHcyfP5+pU6fy6aef8sknn3DLLbcwf/58\nNm3aJOR0hsomtO5O6XK5aAtkn4zNkAhp6NChLF68mMWLFwv/YtCa+NeWexdqysrKuPTSS1m4cCEr\nVqzoFpttlJSU4PP5MJvN3WJ3zZo1zJ49m1dffZXt27cLrUFX2rz11ltZuXIlL730Es3NzUIt7a6y\nGR0dzd/+9jc0Gg2NjY1CAci4uDiqq6upra3FZDIBoNfrsdvtIbUJcOONN3L++ecLn5+MzZAI6cCB\nA3g8HqA1Rai6uhpozQDvqp0BYmJihKqqKpWKhoaGLrfZxurVq5k5cyZRUVHdYtftdgvi0Wq1VFRU\ndLlNq9XKAw88wBNPPEEgEKCurq7Lbbbxv881Li4Os9lMZWUlAE1NTej1+pDZ62gq1Ww2C6sYTmQz\nJGOk5ORkXnzxRaKjoxk3bhwul4sFCxbQ1NTEggULQmHiKO68805efvlldDodEydOpL6+vstttlFc\nXIzZbAZad3zvarszZszgpZdeIjIyUqgu2tU24+LiWLp0KVqtlquvvrpbvtM2pFLpUc/V5XKxePFi\nvvzyS6644ooOa7efCh2V5Jo6dSqLFy9Gq9We0KaY2SAiEgLE8LeISAgQhSQiEgJEIYmIhABRSCIi\nIUAUkohICBCFdJpYLBYefvhhAK688kpuueUWbrnlFm644QY2bNjQ7vMTMWnSpJOy/dVXXwHw5ptv\nsmnTpk6f9+STTzJt2jTq6uq49957ufXWW3E4HADMmzdPmBOcO3cul1566Un51F8Ra3+HgLZ5CL1e\nzwcffACA0+nklltu4Y033uj01iEnu8XIhx9+eNLia7OzePFiKisrueiii0hISODXX38lISGBjIwM\nlEolAAsXLmTq1Kknff3+iCikEPLnKbnGxkY0Gk2Hs+ZtfPbZZ6xZs4aBAwcKuyVs2LCBd955B6lU\nyuzZsxk1ahTXXnstAwYMwGq18thjj+F0OiktLeXVV19FoVDw8ccf8+6776LT6Xj99df57rvv+OCD\nDwgEAjz44INccsklR9lWq9W43W5aWlrQ6/WsWbOGefPmhfah9BPErl0IcTgcQtfuueee4+mnnz7h\nOatXr+azzz7jsccew+PxEAwGWbZsGR9++CHvvfceb7zxBgDV1dW8+OKLrFy5kuXLlzN+/HhSU1N5\n9NFHARgxYgQffPABer2e/fv38+233zJ37lxWrFhxzO1MANLS0mhoaGD37t0YDAbS0tL4+OOPWbRo\nkdDVE+kcYosUQnQ6ndC1a8NisXR4fG1tLfHx8UilUsxmMxEREdhsNsrLy7nzzjsJBoPU19fj8XgY\nOHAgWq0WaM29g/YtYGZmJgBGoxG3283jjz/OP/7xD2pqarj++us79OHxxx8HWsdGs2fP5p133mHy\n5MmsW7eOG2+88YQtqkgropBCSEc/uo7eNxqNWK1WvF4vNpuN+vp6IiMjSU1NZdWqVfh8PlauXIlS\nqaSsrAyXy0Vzc7NQeLOjMVUwGGTt2rXMmzcPpVLJ1KlTueqqqzr0Oz8/n9TUVKGWm0QiEcQq0jlE\nIYWQjn7t6jFSAAAA4klEQVTY27ZtEwbtWVlZPP/88wDI5XLuuOMOpk+fTkpKCgaDAYlEwsyZM7n5\n5ptxuVzccMMNwrGzZ8/GZrPx97//HWhNFl60aBEGg6GdDxKJhMzMTGbMmIFWq2XmzJnH9XvNmjU8\n/fTTQhb9woULeemll457TyLtEZNWewnXXnst69atC8m1nnrqKe644w4yMjJOeOyUKVP44osvQmK3\nLyMGG/opTz/9NPX19cc9Zu7cucI6JJHjI7ZIIiIhQGyRRERCgCgkEZEQIApJRCQEiEISEQkBopBE\nREKAKCQRkRDw/wBnAujYMRd6QAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb66a2b6710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8/2.54,3/2.54))\n",
    "# plt.figure(figsize=(8,3))\n",
    "sns.set(font='Arial',style='white',context='paper',rc={\"xtick.major.size\": 4, \"ytick.major.size\": 0})\n",
    "plt.fill_between(x_scale,(np.mean(rbc_densities[:,rbc_cone_input-697],axis=1)-np.std(rbc_densities[:,rbc_cone_input-697],axis=1)/np.sqrt(len(rbc_cone_input)-1))\\\n",
    "                 /np.sum(np.mean(rbc_densities[:,rbc_cone_input-697],axis=1)),\\\n",
    "                 (np.mean(rbc_densities[:,rbc_cone_input-697],axis=1)+np.std(rbc_densities[:,rbc_cone_input-697],axis=1)/np.sqrt(len(rbc_cone_input)-1))\\\n",
    "                 /np.sum(np.mean(rbc_densities[:,rbc_cone_input-697],axis=1)),facecolor='0.5',alpha=0.3)\n",
    "plt.fill_between(x_scale,(np.mean(rbc_densities[:,rbc_rod_input-697],axis=1)-np.std(rbc_densities[:,rbc_rod_input-697],axis=1)/np.sqrt(len(rbc_rod_input)-1))\\\n",
    "                 /np.sum(np.mean(rbc_densities[:,rbc_rod_input-697],axis=1)),\\\n",
    "                 (np.mean(rbc_densities[:,rbc_rod_input-697],axis=1)+np.std(rbc_densities[:,rbc_rod_input-697],axis=1)/np.sqrt(len(rbc_rod_input)-1))\\\n",
    "                 /np.sum(np.mean(rbc_densities[:,rbc_rod_input-697],axis=1)),facecolor='0.7',alpha=0.3)\n",
    "plt.plot(x_scale,np.sum(rbc_densities[:,rbc_cone_input-697],axis=1)/np.sum(rbc_densities[:,rbc_cone_input-697]),label='rod and cone contacts',c='0.3')\n",
    "plt.plot(x_scale,np.sum(rbc_densities[:,rbc_rod_input-697],axis=1)/np.sum(rbc_densities[:,rbc_rod_input-697]),label='rod contacts only',c='0.6')\n",
    "plt.xlabel('IPL depths [%]')\n",
    "plt.ylabel('Density')\n",
    "plt.yticks([])\n",
    "plt.xlim(50,110)\n",
    "sns.despine()\n",
    "plt.legend(bbox_to_anchor=(0.65, 1.13))\n",
    "# plt.savefig('figures/rbc_density_comparision_on_sac_only.svg',bbox_inches='tight',dpi=300,transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "#### Plot dendritic fields/axon terminals (Supp. figure 3D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Function definition for mosaic plotting\n",
    "def plot_mosaic(selection,layer='OPL'):\n",
    "    im=Image.new('RGBA',(4800,3600),(255,255,255,0))\n",
    "    draw = ImageDraw.Draw(im)\n",
    "    draw.line([(0,0),(0,3599)],fill=(0,0,255,255),width=1)\n",
    "    draw.line([(0,3599),(4799,3599)],fill=(0,0,255,255),width=1)\n",
    "    draw.line([(4799,3599),(4799,0)],fill=(0,0,255,255),width=1)\n",
    "    draw.line([(4799,0),(0,0)],fill=(0,0,255,255),width=1)\n",
    "    del draw\n",
    "    for cell in BC_IDs[selection,0]:\n",
    "        nodes_cell=np.empty((0,3))\n",
    "        for skel in np.where(skeleton_ids==cell)[0]:\n",
    "            nodes=skeletons[skel].item()[list(skeletons[skel].dtype.names).index('nodes')][:,:3]/[16.5,16.5,25]\n",
    "            if nodes.shape[0]<2:\n",
    "                continue\n",
    "            nodes=np.dot(M,nodes.T).T\n",
    "            edges=skeletons[skel].item()[list(skeletons[skel].dtype.names).index('edges')]\n",
    "            soma_pos_cell=soma_positions[BC_IDs[:,0]==cell,0]*16.5\n",
    "            draw = ImageDraw.Draw(im)\n",
    "            if layer=='OPL':\n",
    "                nodes=(nodes*[16.5,16.5,25]+[0,3250,15000]).astype(int)\n",
    "                for i in range(edges.shape[0]):\n",
    "                    if (nodes[edges[i,0]-1,0]<soma_pos_cell) and (nodes[edges[i,1]-1,0]<soma_pos_cell):\n",
    "                        draw.line([tuple(nodes[edges[i,0]-1,1:3]/[25,25]),tuple(nodes[edges[i,1]-1,1:3]/[25,25])],fill=(0,0,0,150),width=5)\n",
    "                del draw\n",
    "                nodes=nodes[np.unique(edges).astype(int)-1,:]\n",
    "                nodes_cell=np.concatenate((nodes_cell,nodes[nodes[:,0]<soma_pos_cell,0:3]),axis=0)\n",
    "            elif layer=='IPL':\n",
    "                nodes=(nodes*[16.5,16.5,25]+[0,5000,22750]).astype(int)\n",
    "                for i in range(edges.shape[0]):\n",
    "                    if (nodes[edges[i,0]-1,0]>soma_pos_cell) and (nodes[edges[i,1]-1,0]>soma_pos_cell):\n",
    "                        draw.line([tuple(nodes[edges[i,0]-1,1:3]/[25,25]),tuple(nodes[edges[i,1]-1,1:3]/[25,25])],fill=(0,0,0,150),width=5)\n",
    "                del draw\n",
    "                nodes=nodes[np.unique(edges).astype(int)-1,:]\n",
    "                nodes_cell=np.concatenate((nodes_cell,nodes[nodes[:,0]>soma_pos_cell,0:3]),axis=0)\n",
    "            else:\n",
    "                print(\"Layer has to be 'IPL' or 'OPL'\")\n",
    "                return im\n",
    "        if nodes_cell.shape[0]>2:\n",
    "            nodes_cell=nodes_cell[:,1:]/[25,25]\n",
    "            hull=ConvexHull(nodes_cell)\n",
    "            draw = ImageDraw.Draw(im)\n",
    "            for simplex in hull.simplices:\n",
    "                draw.line([tuple(nodes_cell[simplex[0],:]),tuple(nodes_cell[simplex[1],:])],fill=(0,0,255,255),width=5)\n",
    "            del draw\n",
    "    return im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_mosaic(rbc_rod_input-390,'OPL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plot_mosaic(rbc_rod_input-390,'IPL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plot_mosaic(rbc_cone_input-390,'OPL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plot_mosaic(rbc_cone_input-390,'IPL')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}