CBC5_classification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CBC5 classification\n",
    "This notebook contains the reclassification of the subtypes of CBC5 based on the density profiles of their axon terminal systems in the IPL. The resulting classification, together with the CBC1/2 classification from the corrigendum of Helmstaedter et al. (2013) is written to a file to be used by all other notebooks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat\n",
    "from scipy.spatial import ConvexHull\n",
    "from PIL import Image\n",
    "from PIL import ImageDraw\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from sklearn.mixture import GMM\n",
    "from sklearn.decomposition import PCA\n",
    "from shapely.geometry import Polygon\n",
    "from descartes.patch import PolygonPatch\n",
    "\n",
    "%matplotlib inline\n",
    "matplotlib.rcParams.update({'font.size': 14})\n",
    "sns.set_style(\"whitegrid\")\n",
    "matplotlib.rcParams['font.family'] = 'sans-serif'\n",
    "matplotlib.rcParams['font.sans-serif'] = ['Arial']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "BC_IDs=np.loadtxt('data/cell_IDs.csv',dtype=int,delimiter=',')\n",
    "BC_IDs=BC_IDs[(BC_IDs[:,3]>57)&(BC_IDs[:,3]<72)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#CBC1/2 classification from corrigendum\n",
    "CBC1=[393,400,403,410,412,414,415,417,418,420,421,422,424,427,430,431,432,433,436,438,440,441,444,445,447,449]\n",
    "CBC2=[390,391,392,394,395,396,397,398,399,401,402,404,405,406,407,408,409,411,413,416,419,423,425,426,428,429,434,435,437,439,442,443,446,448]\n",
    "BC_IDs=np.concatenate((BC_IDs,BC_IDs[:,3].reshape(-1,1)),axis=1)\n",
    "BC_IDs[np.in1d(BC_IDs[:,0],CBC1),4]=58\n",
    "BC_IDs[np.in1d(BC_IDs[:,0],CBC2),4]=59"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "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": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "volume_density=np.loadtxt('data/density_data_BC_flattened.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "start_ipl=[]\n",
    "for i in range(volume_density.shape[1]):\n",
    "    try:\n",
    "        start_ipl.append(np.min(np.where(volume_density[:,i]>0)))\n",
    "    except ValueError:\n",
    "        continue\n",
    "end_ipl=[]\n",
    "for i in range(volume_density.shape[1]):\n",
    "    try:\n",
    "        end_ipl.append(np.max(np.where(volume_density[:,i]>0)))\n",
    "    except ValueError:\n",
    "        continue\n",
    "volume_density=volume_density[np.min(start_ipl):np.max(end_ipl),:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cbc5_volume_density=volume_density[:,(BC_IDs[:,3]>62)&(BC_IDs[:,3]<66)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cbc5_volume_density2=np.copy(cbc5_volume_density)\n",
    "cbc5_volume_density2[:1500,:]=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "on_sac_volume_density=np.loadtxt('data/density_data_ON_SAC_flattened.gz')\n",
    "off_sac_volume_density=np.loadtxt('data/density_data_OFF_SAC_flattened.gz')\n",
    "on_sac_volume_density=on_sac_volume_density[np.min(start_ipl):np.max(end_ipl),:]\n",
    "off_sac_volume_density=off_sac_volume_density[np.min(start_ipl):np.max(end_ipl),:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Type 5 subgroups by density profiles (according to Greene, Kim et al. 2016)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "density_PCA=PCA(n_components=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PCA(copy=True, n_components=3, whiten=False)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "density_PCA.fit(cbc5_volume_density2.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot PCA components (Supp. figure 7A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI4AAAB/CAYAAAA5M/q0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl4lNXZh+9ZM5ksk8m+r4QQyE4ghB2sIq1aUFSUxYq2\ntraVSmldKq0Lat2qlWpVUKwoILYV5WtlVWQNZAOyQIDs22Sy7zOZ7fsjEIghycwQk0k693VxXTrv\nOed93slvzjnvOc/zHIHJZDJhx46FCEfaADujE7tw7FiFXTh2rMIuHDtWYReOHasQW1ohMzPz+7DD\njg0wefJks8taLBxLbzCcnD17lujo6JE2o19s2T5LOwT7UGXHKuzCsWMVo0Y4F8ubKKxoGmkz7FzC\nqjnOcHMou4JXPu4eg5/9WSqJUd4jbJEdm+9x1A0dvPZJJqmxfsRGePL0pjRe2ZJBY6tmpE37n8bm\nhXP4VCUikZDfLE3kmZ9N44554zh1oZa1bx6mVNUy5Pfr6Ooc8jbHIrYvnNOVJEf7IJdJkIhFrPzh\nRF5bPRuJSMDr27IYys39A4VH+Mnna8isyhmyNscqNi2cqro2CiuamRUf0OtzXw8nHvxxLIUVzRRW\nNA/Jvbr0XXyauwuAlw6/zbEy+0LnQNi0cE6dr0UogMnRfSfDieO9UDhL+SarfNB2VG21lDdXDVjm\nXF0hTZoWfj5lBV5yd944volH//sMTZ1DI8yxhk0L51xJA6H+CuQySZ9rIpGQ2YmBHMquxGAw9ttG\np07Dmq+e5be7n+O99E8wGq9d9nx9EXKJI3PDprHhludIDoinslXFL/7vD+iNhiF7prGCbQuntJEJ\nIcp+r89NCqSpVUv2+dp+y3xTfAy9UY/CwYX9RUdY+tkvOVh8vE+5grpCIj3CEAqECAVCfj/z59yX\nsASD0cDZ2gtD8jxjCZsVTnObluq6diaEuvdbJjLIjQAvZ77J7H+4OlWdR6xPFBsXvcxjsx7GQ67k\nnfSPqWpR9ZQxGA2crysmyjOiV92F4+fh4ajkZMWp63+gMYbNCqegtBGACSH9C0cgEHDDlCCO51TT\n2tHV57rRaKSgvohor/EATPaP5a8Ln0YiknC8PKunXGZVDp16DYl+k3rVFwqETAmM52TlKYym/ofD\nwTAZDHQ1NFpd3xaxWeGcK21A4SzF10M+YLkbpgRjNJr4OqNvr1PVVkOnTsN4j7Cez6RiKYm+kzhR\nkQ2A0WRkV8F+xrmHEuEe0qeNKQHxNHY2U9ZUCYBBq0VbW4fJZOr5NxgFr/yF9PsfRFdcMmjZ0YLN\nbjmcuVDHxDAPBALBgOXcXWWkxPiy+3gJkbf69bpW0lgBQKgyqNfnM0Om8OrRdymoKySrKpeCukKe\nmP3La7Z/WUxFjeUEiNzI+uUj6Fu6Fx4FIhGK2Bii1z2JUHztr7L228PUH08DoPXdTWR+8X/4L7oN\nv4ULBv4CbByb7HGa27ScL28kOdrHrPI3TwulQt1GUXXvVd+SpnI8HJW4Ojj3+jzZP44AV1/+eOA1\nPj+7m7tibiXRL+aabcsljvg6e1HSWE5Dejr6lhYcA/wBcAwMoOnUaUr/seWadU1GI2Vbt6OIjSH1\ns23If7QQQ2cnJR98SFfj6B66bFI4WQVqTCaYPMG8zcz4SC+CfV3YnVHXa+gobiwnRBnYp7xQKOTR\n1AeZFpTEz6es4I6JCwdsP1QZRHFTOTV79uM6MZqktzcw44t/kfjm6wTetQTVnn0YtNo+9ZpzctGo\nVATdfSdCqRTZrBkkbngdhEJUu/ea9Wy2ik0KJyO/hnB/BR4KR7PKC4UC7r1pAqVqDar6DroMOnae\n3UOe+jyx3lHXrBPsFsCj0x9kfvj0QYfDcGUwNdWltBYU4P2D+b2uec6YjlGrpTknt0891Z69yPz9\ncY25MumWKBR4TJtK3eEjQ7pdMtzYnHAMBiNZBeprrhYPRFykJwC5hXVsOfUvtp7Zib+LDz+ImHXd\nNo33CMOvvBUAZVJir2vykGDkoSGcfe4FVHv3UfTeJhozs+hqaqIh7SS+C27sI0zPGdPprKyio7T0\num0bKWxuclxQ1khbp44p0b4W1XORS/FWSDlf3kCu9AwSoZhn5q/BQSy9bpsi3EMJq9JhCPRGquy9\nICkQCAhb9RMuvPkWhW+9A0D1f75CIBaDQID3/Ll92nNLTEAkl1P77WGcQkOv276RwOaEk3G2Bhe5\nhPEDrBj3h5+HA2frCqhXNrL+ht/h8p1JsbU4iKUE1RlQx127Pbf4OJLfe5u2wiIcvL0p3fIxjekZ\nBC+/F4mra5/yQokEj2lTqfz3TtynJOM60TYd2AfC5oSTnl9DYpQ3IuHA845r4eshpqDjFAEuPkRe\ntXZzvRi0WmQdOkolHf2WEYhEuIyPBCDy19d+tb+akPtW0l5aTs4TTxH56CN4z50zZPYOBzY1xymu\naqakuoXpsf5W1Tc6qxE4NXNrxK2DTngtQVvbvRdWKdFQ3zE0r9FSNwWTnl6H8/hICv/+HvqO0eVA\nZlPC+ep4CW7ODkydZNn85jL1gjKMHc7INH6DF7YAbY0agBZnEbk1BUPWrsTVhai1azBqtdQdOTpk\n7Q4HNiMcnd7AwcwKbpgShERsuVl6g56ijlIk7f4UVg5tNISmRg1CIW5+geSozw1p2zIfb5TJSVT+\n698Ydbohbfv7xGaEk1/cQKdWz8zvePuZS476HFpjF+HOUeQXNwypbVq1GgdPT+L9J5FdnYdhiP1z\nQpYvQ6OupWbfgSFt9/vEZoSTXaBG4SwlPEBhVf208myUUleSQyM5X9ZIl27o/riaGjUO3l5MC0qi\nVdtGnvr8kLUN4BQagjIpEfXXB4e03e8TmxFOVoGahEhvhFa8TXXoOkkrzyLGbTyxEV7o9EYKyoZu\nL0hTU4PM14dx7qF4yt1Ju7SzPpR4zppJ24ULdFarBi9sA9iEcBpbNBRXtZAY5WVV/czKbn+aye4x\nhAUokMvE5BbWD4ltJqORzvIK5IGBCAQCpgUmcrIie8iHK4+UKQilUuoOHxnSdr8vbEI4l10/rY3Q\nPF6RRYR7CEqpKyKhgIlhHuQW1g2JbZoaNcauLuTB3a4Z04KSaPkehiuRoyPuU6dQe+jwqNjDsg3h\nFKgJ9XPF3VVmcd1OnYbT1XmkBiX1fBYb4cG5kgZ0euu89kwmE+U7/knen56l5IMPQSjEKSIcgEiP\nMLydPDhSlm5V2wPhOXsWneUVo2IPa8SFozcYyTxXY7YLxXfJqs5BZ9STEnhl83FSuAddeiMXy617\nLW/Kyqbsk200nTpNw8l0/G/5IVI3N6B7b2p6cDInK06hM1z/63NDi4aDWRW0d+pQJiUgcnKi9pDt\nD1cjLpxT52tp7dAxM8G61/DMqlxCFAH4OF+ZH0UEuuEgFZFbZN1wVf2f/+IUFkrihtcJf+inBC+/\nt9f1GcHJdOg6Oa3Kt6p96F63+mBXHj//8wFe+ySTZ99PwyAQ4ZGaQv3R4zY/XI24cA5lV+Dn4USE\nFa/hRpORHNVZ4v0m9vpcLBISHeJObpHlE2SDRkPTqTN43zAPeXAwfj+8GZGDQ68ywYoAfJ29yKrq\n64NjLv/+5iKfH7zIjDh/nrhvCufLGnl5SzruKVPRqFR0lg8eaDiSjKhwquvaOZRdyYJpIVbtLV2o\nL6ZZ20qC78Q+12IiPDhb3DBgsN61aC8pxWQw4DJhQr9lBAIBMT4TyFVbt/1QoW5l+74CFs2JYPXS\nRKbH+bN2eTJpuSrSNW4IHRyoTztpVdvDxYgK5+uMcqQSET+aYd1O9q6C/Xg7eTDxUvjL1cREeNKp\n1VNUZVkIr6a6GqDHr7g/YryjULXVUtdu2Sq13mDk3c9zULrKWL7wijvFjDh/5icHsfG/BcgmxdBw\nwi6ca2IymTh8qoKpE32ROVju3dGh6yS7Kpebxs1BKOz7GJFBbkjEQovXczTVKiQKBWL5wGE5Md7d\nYj2lyjO77ezCFn76wn5yLtbxi9vjcJCIel1/aHEsrnIJGSYf2i4Woq0fmrWo74MRE05OYR2Vte3c\nmBJsVf1DJSfQmwxMD752BlSpRERUiNJi4XRWq5D5Db477ypzIc4nmgNF197VNhhNnMxTcexMFVqd\ngUPZFWz7RkWIrwuvrp7NlIl97yGXSbj/1knsbXLBJBDQcDLDItuHkxFz5Np9vBR/TyfixnlaXFdv\n0LPr3D5SAhPxlPcf6RkT7smuI0UYDEZEIvN+I5pqFY6B5r3h3ThuFq8dfY+L9SWM8wjtde3df5/h\nq+MlADg7SujQ6okLc+ZPD04bcD43OzGQg1nBqOt9cTueZrPxVyPS4zS1ajmeU8WCaaFWTYq/Lj5G\nXUcjt0cPHNaSHO1Ne6eO0xfMfy3XqKpxNKPHAZjsH0eAiy8vHX6bjRlbe9Z1aho62J1Wwr0LJvDn\nX87EXSHjtlnh3DPPz6znvTk1lAxZGM2nz9BeWma27cPJiAhnf3oZ0B33bSkavZbP83eTEphI6DVi\npq5mfLCSAC8nDqSb9+XrWlvRt7aZNVQBiIUiHkldhVziyL7Cw3yW9x8A9p0oRSYVsWhOBJPCPXjr\nd/N54LYYs91h48Z5ctEtAoOTCxc3vEXdMdtb1xl24RiNJvaklTA9zg+Fs8PgFa5Co9ey9fROmrWt\n3Bv340HLCwQC5iUHkZZbbdb2g+bSzrTMz3wPwjBlEG/88GnunPQjvjy3j+LGCvanlzE7MRBHKyb9\nADKpmInjfUiPXoCxq4uCl14lb93TFL23ia4m20j0NOzCOX2hFlV9BwtTQy2qZzKZeOXI39l98SC3\nRN2Ar4t5WxRxEV506Y0Um/Faflk45g5VlxEIBCyKXoCPkyfvpm2nvlnDvMlBPXZbQ3K0D4dbnIh8\n6WWifr8Wg0aLavdeCl56BUPnyPsnD7twdqeVEOTjzKRwD4vq5dScI6emgCS/GJZM+pHZ9SICFYhF\nQrO8AjurqxG7uCB2tjysRiKScFfsLRS1FOLk3saEUHeOlWWy6vPfsvJfv+GNY5vQGPqGCffH1Em+\nGE0mjp6uxHNGKvGv/pmY9c/QVlTM2RdewtjVN63LcDKswmlo0XAiV2XxpNhkMrHtzBeMcw/lsVkP\nIxX1Te3WH1KJiIgABRfKB3fs0pj5Kt4f0wKTkOhdcAy7wKe5X/DG8U1EeoQR4R7CsfJM3jj3IRUt\n1Wa15a2Ukzjemz1pV3bKXSdGE/3kY7Tkn+X02sfQ1o3cOs+wCmffyVKEQgHzky2bFGdX51HYWMrS\n2NusegsLD1SYlZ1UU63C0YL5zXcxGKCzJJJ2UQ07z+5holckv0l9kD/Ne5T1N/wOsUDMP/P+a3Z7\nN6eGcKG8iZyrfIvc4uOIf/UldC0tFL+/2Wpbr5dhE05zm5Yvvi1kblIgLnLLwnIPlaQR6OpHrE//\n+0cDERHgRlVdGx2agd0gNKrq6+pxLpQ3oWvwYm7APNZM/ylPz1+DXNqdOGG8ZzhzfVI4VpZB2lXZ\nwAYiZZIf4f4Ktu3pvSfmFBZK8NK7qT92nPaSEqvtvR6GTTgf/fcsBqOJFQstC3dt1baRXnma2aEp\nVgfZRQQqMJmguKr/TOy6lhZ0zS04+lsXDAiQX1yPg1TMz1KXMO0qx7LLTHafxNTABN5J/5gO3eAT\nXKFQwOK5EeQU1lH2nSzy3jfMw8Hbi+IP/jEiYTXDIpyC0gb2nSxl2YIJKC308vu25ARGk5G5YalW\n3z/E1xWxSDBgvFV7SfdcQh7aN52bueQXNzAhRIm4n1VqgUDATxLupEPXSUblGbPanBHvj6uTtGcV\n+jJCiYTwnz5AS14+Fze8PezrPN+7cDo0Ov76aTahfq4W7YK3dbVzsb6Ef+d/xZSABNxkfYP3zUUi\nFhLi5zrgPKe9uBihVIrczO2G76LTG8gvrmdS2MBvi55O7oS6BZJdbZ4vj0Qs4sapwRxIL+8z1LpP\nnULkI7+i9ttDVOz4p1V2W8v3Kpy2ji6e3phGbWMna5dNNnu/qOj9Dzhx3wP89dP1SE1CVkTdct22\nRAS4UVTZv3DaLhYhDwlBIBL1W2YgMs+p6dDoSY0bfKhL9IvhlCrfbNfTH80IR6sz8MWhIqDbn+dA\nehkHsyqQJE0l8M47KNu+g+Y86z0SLcWqpc1Psz7HrUWPUifCVSzHVeyIIjAEkVyO2MmJRp2AvSdL\n2ZN1li6NgD89OJNgX1dMBgMNGZkIpVIcvDwp2/opLbl5+C+6jYAf34pAJKIpP5/qL/+DDLhrXxNC\nUxMXP12Ly5uvI/Ox/pyqiEAF+9PL0OoMfdwZANqLilDETLpGTfM4lF1JsK8LoX6D94yzQqby+dnd\nHC/PYnZoyqDlvZSOLJ4Twfa95zh8qpLymtaea0IBxIcH8wP/YApee4OInz0AdCeAEkqvPzdQf1gl\nnMBnPgZAB9Rf+nc1eqEAsYeYWc4idGIBuz7eSbZfOBOy1XRVXQk4kygUuMZMpPSjj6na+SUiJzkd\nDfWoPcT4//4RfL/NR1tTQ2NmNpk/+wX+P76VsFU/sepBIwIUGI0mSqqaifpO7mSDRkNnZRX+t1nX\ns3VodJzIU3HXDyLNKh+o8CPWJ4rdFw6aJRyAFQujcXQQU6Zq5Z4bo0iI6g48PJGn4lB2Be+Jk7i3\neQ+6F18GQOjoiFSpxHVCFCJ3d9rbOtFWVCL198dv5jQY58fJilMcLcugpr2Oh4OWWvTMVgknLXI+\n9cioEgjQSbowOXTiiRqZyICLSUCASERIYwfBXUaMrVpMRQ0Iss5Q7uWA52+WkRSaSGdZBe5TkxE5\nOtJWWETtocPUVZVSKmhAeu8iUifMgAkzANCoVJRt20HVF7twS0xAmZhgsc2h/grEIgEFZY19hNNZ\nVQUmE/Jg63yD0nJVdOkMzEkceNP1ahZGzuPlI+9wvq6I8Z7hg5YXiYTcfWPffIYLU0NZmBpKffNk\n9qclcCzzIpVVjUS1lyFv7yT02EmkBh0mgYBGR0c8C/Jo3fsVWomAdjcxwsRIJPJwsHC/2Srh/PbV\nK4mD9AYjeoMRiVjU7+6vrrmZhkY1JyoPsr1iH5MoxVEso/noceaEpnBjxGw6fFx5Zu8LBE9LYl3q\nkl71Zb6+RP7m13SUV1C+bQduCfEWv5o7SEREBnU7dt02q3fq/c5K89xF++ObzHImhrnj6+Fkdp0k\nv1j8XLzZemYnf5r36HXn8yloziPdtBuPmW5MMElRNzvS2NXBEcMVm4QIcZZ4MqnOFa9qIx4XK1i4\nN58uaTHETbXoftftyCUWCft9/byMRKHAR6FgTWgkR8vS+U/B14iFIpylcjZlbud4eRZVLTU4SeX8\nftYvkFxjS0EgEBCy/B7yn1mP6qs9+C5cYPGXHTvOk6+OFWM0mnrFqGuqqhA7OyN2cbGoPeg+U+v0\nhVoeviPeonpCoZCfJN7Ji4fe4nh5Vr+ejIOhNxr4NOdLvji3l0ne4zGaTHToW/BWOhHlGMCskClE\nuIdiMBqQiR16fWdGnY7642loVDWoLbzvsHsAzgiewozgKT3/f6wsg/ezPsVb7sHamQ/hLO3/V+uW\nmIDnzBkUvbuRxoxMov/wuEVvQXERnuzYf55SVQth/lfCcTqrqpD5m+dkdTXNbVqe3XQCN2cH5k42\nf5i6TKJfDPG+0Xye/xWpQUkW319n0PHCob+RX3uBpbG3sSh6AUJBPz/ia/wYhRIJXrO7s7KqLTyw\nfsRzAE4PTiY1qPvXNtgXJxAIGL/2UdwSE7i44S1Ue/fht/Bms+8VFapELBKQc7Gut3Aqq3AMsHz9\n5oNdedQ3d/LcQ9ORSa37KhdGzuPPh9+msKG0j/vpYGzP+ZJztRdZN+cRYqzcjrGWEQ/Ig25BmPtr\nEwgE+PxgPl7z5lL2yTZ0ra2DV7qETCpmfLCy16ahyWSio7wCeZBlPUZLexeHsitZemPUgEcjDUa8\n70TcHd3YV3jYononKrLZVbCfu2NvG3bRgI0IxxpCVy7HqO2iZu9+i+rFjfMi52Jdj0egtrYWo0bT\nk43CXA5mlmMymZhvhfvr1YiEIhZGzuPbkjRUrWr0RgOqtlry1ReoaK5Gq++iskXVc+xRfUcjrxx5\nh9eOvsdk/1hui7rxuu5vLSM+VFmL1F2JW1IidUePE3jHYrPrpUzyZfu+AnIL60iM8qajrDvU1lLh\nHMgoJznaB6WL5Rk2vsvNkXPZVbCPvxzbSLuuk9r2vn42XnJ34n0ncrQsA6lIwsNTVzIrZOo1Y8qG\ng1ErHADPGamcf+2N7oxZPuadNBMRqMBTISMtt7pHOEIHBxy8zV+VbmzVUFTZzI9nRwxe2AwcxFJW\nJd3NexlbEQqErJ3xED7OnlS3qqnvaMRDruRk5WmOlWcyPTiZZfGLBnyJGA5GtXCUyckIJBLqj6UR\nsHhw53XoniNNi/HjeG41Dy2Oo6OsHMfAQAQW/HIvB/nFRlgeE9Yf04OTSbh0Qp9c0u3DE+J2Zd51\nLTeNkWTUznEAxHJHlEkJ1B3tezjrQEyL8aO+WcPFiiY6yspxCrFsmMo6p8bP0wkvpXmn25iLXOLY\nIxpbZ1QLB8BjeiptFy6gUZu/hDUpwgMnRwlpOZV0VlTgGGS+cLp0Bo7lVDEz3nqHr7HAqBeO+5Rk\nBGJxz/GF5iAWCZke68fJo/kYtVqLJsbp+TV0aPTMTbJ8wW8sMeqFI3Zywi0hnvrjJyyqt3juOKQN\n3b2UJcL5JrOc8AAFwb7WO5aNBUa9cADcU6bSeq4AY1ub2XWCfFxIVBrQCSWI3c2L8Tp9oZYTeSpu\nSrHevXSsMDaEMzUZAN05yzJkRTlqqZUoOFtqRsyVVs/fPjvFxDB3i6NQxyJjQjhSNzdkfn7oK6ss\nqidpqKHF2Z2jpwev9/Huc9Q3a/j1XQlWZX8fa4wJ4QA4hYVgqDIvShLAZDCgqazCfVwYBzPL6dTq\n+y177EwVXx4u5J6bogj0ttz1YiwydoQTGopeVYPJaF6yyM5qFcauLuJnxNGp1bPvZN+k1FqdgXc/\nP8NLWzJIiPRi8dxxQ232qGVUrxxfjVN4GGi1aNVqZL6DR2O2F5cAEBg/gZllYj4/WMjN00KRXnJk\n79DoePytI1So21ixMJrFcyLMjtL4X2DMfBNOYd0xW21FxWaVby8uRqJUInVz4+4bx9PYouGzAxeA\n7iOs139wkqq6dl761UyWzI+0i+Y7jJkeR+quRODsTHthEZ7TB4/6bMk/i/Ol8xmCfV25fd44duwv\noEOr49T5WspUrTz902lEBll+GvH/AmPmZyQQCBCHBNN06nS/ZToqKlEfPERndTWtZ8/hOXN6z7V7\nbopifnIwXx4qQqcz8ssl8UyeYN6O+/8iY6bHAZDGxdC2bcc1PfpMRiO5f/gjuqYr8eMeqdN6/lsi\nFrF6aSLLF05A6SKzv3IPwpjpcQCkE6MROTmh/uZgn2tNp06ja2rCJao7sXXIyuWIZH2dsDwUjnbR\nmMGY6nEEEgke01KoP3qMkBXLevkxqw98gzw4iNiXXhjSM8n/VxlTPQ6A58zpaFQ1tJ2/0POZyWSi\nOTcPZfJku2iGiDEnHLf4OBy8vSjbvqMnZ0xHSSm6pqbrSipgpzdjTjgCkYiwVffTlJXNuRdeoj7t\nBMUffIhE6YYiLnakzRszjKk5zmU8UlOIfHQ1Re9touFkOgiFTPzjHxBKzM9WamdgxqRwALznzsZr\n9ky6GrtdJhw8LMurbGdgxqxwAARCoV0w3xNjbo5jZ3gQmCxMV5lpYVYDO6OHyZPNT7VisXDs2AH7\nUGXHSuzCsWMVduHYsQq7cOxYhV04dqzCrAXAtWvXMn/+fKqrq6msrKStrY0nnngCpXJk3Sqrqqp4\n+OGHiY6OxtvbG1dXV5uyr7KykrfffhtnZ2dcXV2RyWQ2Y98nn3xCbm4uOp2OrKwsli1bZpFtg/Y4\nH374IU5O3Ul8MjIy+OMf/8gdd9zBjh07huYJroP09HS8vLwASEhIsDn7Nm/eTHBwMC0tLcTGxtqU\nfcuWLePFF1/Ez8+Pv/zlLxbbNmCP8/XXX+Pi4kJCQgJGoxGPS8v3vr6+1NbWDs0TXAdxcXHMmDED\nT09P7rvvPoIupSuxFftKS0tZsmQJDz74IPfffz/BlzK324p9RUVF6PV6/P39Lf7bDiicXbt2oVAo\nKCrqPrXkcs+jUqnwtiD12ffF2bNnSUjoTs8vk8lQX8qRYyv2eXl54eTkhEgkQiaT0XTJ39lW7Nu6\ndSurVq3Cw8PDYtvMWjneuXMnUqmU+vp6ioqKaGlp4ZlnnsHZitNyh5K8vDw2btyIp6cnUVFRaDQa\nm7KvsLCQDRs24OLiQkpKCo2NjTZl3/3338/mzd3nem7ZssUi2+xbDnaswv46bscq7MKxYxV24dix\nCrtw7FiFXTh2rGLUCKeyspJHHnkEgJtuuomVK1eycuVKli5dyoEDB3pdH4xFixZZdO+dO3cC8Le/\n/Y1vv/3W7HqPP/44d955J/X19Tz00EPcd999tF1KcLlu3Tq6uroAeOqpp5gzZ45FNo00o8pZ/XIU\npqurKx999BEA7e3trFy5kg0bNlh0dJElbNmyxWKxXb7Piy++SHV1NTNnziQgIIBjx44REBBAZGQk\n0kun9K5fv54lS5YM0pptMaqEc5mrl56am5uRy+UMthz12WefsX37dsLCwtDru/P9HThwgI0bNyIU\nClm9ejUpKSnceuutBAUFoVKpWLNmDe3t7RQXF/P6668jkUjYtm0b77//Pi4uLrz55pvs2bOHjz76\nCKPRyK9+9Stmz57d596Ojo5otVo6OztxdXVl+/btrFu3bmi/lGFm1AxVV9PW1tYzVD377LM8+eST\ng9bZunUrn332GWvWrKGrqwuTycTf//53tmzZwgcffMCGDRsAUKvVvPzyy2zevJl33nmHBQsWEB4e\nzqOPPgpAUlISH330Ea6uruTl5bF7926eeuopNm3a1CPI7xIREUFTUxPZ2dkoFAoiIiLYtm0bzz//\nfM/QNdoYlT2Oi4tLz1B1mcrKyn7L19XV4evri1AoxN/fH6VSSUNDAxUVFTzwwAOYTCYaGxvp6uoi\nLCysZ7mQiX4iAAABeklEQVRdq9UCvXu4qKjuI5zd3NzQarWsXbuWt956i9raWu66665+bVi7di3Q\nPbdZvXo1GzduZPHixezatYt77rln0B7T1hiVwunvS+7vczc3N1QqFTqdjoaGBhobG3F3dyc8PJx/\n/OMf6PV6Nm/ejFQqpbS0FI1GQ0dHB6JLB8X2NycymUx8+eWXrFu3DqlUypIlS7j55v7PCM3Pzyc8\nPByxWNzT7mVxjjZGpXD6+0OeOHGiZ5IZHR3Nc889B4BYLGbVqlXcfffdhIaGolAoEAgErFixgmXL\nlqHRaFi6dGlP2dWrV9PQ0MBjjz0GQEhICM8//zwKxZUDYC+fIxoVFcXy5ctxdnZmxYoVA9q9fft2\nnnzyyZ6d8vXr1/PKK68M+Ew2i8lOL2655ZYha+vxxx83nT9/3qyyt99++5DddzgYlZPj0cSTTz5J\nY+PAZ0U89dRT1Nf3PYfTlrG7VdixCnuPY8cq7MKxYxV24dixCrtw7FiFXTh2rOL/ASWwrvoxQaT5\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbc8b433400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5.4/2.54,4/2.54))\n",
    "sns.set(font='Arial',style='whitegrid',context='paper')\n",
    "plt.plot(np.linspace(-17.423,119.651,3576),density_PCA.components_.T)\n",
    "plt.xlim(40,70)\n",
    "plt.xticks([40,50,60,70])\n",
    "plt.xlabel('IPL depths [%]')\n",
    "plt.yticks([])\n",
    "# plt.savefig('figures/type5_PCA.svg',bbox_inches='tight',dpi=300,transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cbc5_params=density_PCA.transform(cbc5_volume_density2.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cbc5_GMM=GMM(n_components=3,random_state=198734)\n",
    "cbc5_GMM.fit(cbc5_params)\n",
    "cbc5_prediction=cbc5_GMM.predict(cbc5_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cbc5_prediction=cbc5_GMM.predict(cbc5_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAESCAYAAAC2KnFgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHkhJREFUeJzt3WtsVNfZ9vHLnB4g4yP2jI1eiTaQqRyimEcQocCrElVC\nJCEoeaWGVqEQq23iKEpgpFQEW3EoaoVFVVq5jRILYmErUptYqTAvTTEqisoHxyXBAXIaySTocYqx\nPTZkfIoHE9jPB2NiE4y957TX7Pn/PmXbm/G9R+1cs9a+19oZlmVZAgDAEDOcLgAAgPEIJgCAUQgm\nAIBRCCYAgFEIJgCAUQgmAIBRZtn9B+3t7aqurlZeXp6WLl2qS5cuqaOjQ4ODgyovL9fIyIj27Nmj\nnJwcLV68WJs2bVJtba3tc3JzcxNxvQAAw2XYXcf06aefKicnRz6fT2VlZZozZ45ee+01vf/++zp1\n6pQuX76sH/7wh1q2bJnKysr0hz/8QS+88IJqamqmdc6JEyd0+vRplZWVJeqaAQAGsz2Vt3TpUs2a\nNUtlZWVauXKlFixYIEny+XwKhULq7e1VUVGRJCkrK0sDAwPKz8+f9jmFhYXq6emJy8UBAFKP7am8\nzz77TAsXLlRtba22bt2qq1evSpK6urrk8/l07do1dXZ2yufzqa+vT16vV+Fw2NY5Xq930r/f2toa\nzXUCABJo+fLlcXst21N5H330kV5//XUVFhZq9uzZKiws1Llz59Tf369du3YpEomoqqpKHo9H99xz\njx5//HG98cYbts/xeDy3/Putra1xfQPiLRgMqri42OkyJmVyfSbXJlFfrEyuz+TaJPPri/fnsu0R\n07333qs//elPk/7e4/Fo7969E362efNm2+cAANIT7eIAAKMQTAAAoxBMAACjEEwAAKMQTAAAoxBM\nAACjEEwAAKMQTAAAoxBMAACjEEwAAKMQTAAAoxBMAACjEEwAAKMQTAAAoxBMAACjEEwAAKMQTAAA\noxBMAACjEEwAAKMQTAAAoxBMAACjEEwAAKMQTAAAoxBMAACjEEwAAKMQTAAAoxBMAACjEEy3EQ4P\nqLS0SatWHVdpaZP6+gacLgkAXG+W0wWYLBBoVn39OkkZammxJB1VXd2DTpcFAK7GiOk22trmScq4\nfpRx/RgAkEgE0234/cOSrOtH1vVjAEAiMZV3G9XVqyUdVVvbPPn9w9ePAQCJRDDdRnZ2JveUACDJ\nmMoDABiFYAIAGMX2VN6pU6f05ptvyuPxKC8vT3PnzlVHR4cGBwdVXl6ukZER7dmzRzk5OVq8eLE2\nbdqk2tpa2+fk5uYm4noBAIazHUz9/f3auXOn5s+fr1/84heaM2eOXnvtNb3//vtqaGjQ5cuXtWXL\nFi1btkxlZWV67LHH9MEHH6impmZa55w4cUINDQ0qKytLxPUCAAxnO5jWrFkjSaqpqdGGDRt08uRJ\nSZLP51MoFNKVK1dUVFQkScrKytLAwIDy8/OnfU5hYaF6enpivzIAQEqyHUxDQ0PavXu3NmzYoPvu\nu0/Hjh2TJHV1dcnn8+natWvq7OyUz+dTX1+fvF6vwuGwrXO8Xu9tawgGg3bLTppIJEJ9UTK5Non6\nYmVyfSbXJplfX7xlWJZlTX3atyoqKvTll19q4cKFmjlzpu6++26dO3dO/f392rVrlyKRiKqqquTx\neHTPPffo8ccf1xtvvGH7HI/Hc8u/39raquXLl8fl4hMhGAyquLjY6TImZXJ9JtcmUV+sTK7P5Nok\n8+uL9+ey7RHT7t27b/t7j8ejvXv3TvjZ5s2bbZ8DAEhPtIsDAIxCMAEAjEIwAQCMQjABAIxCMAEA\njEIwAQCMQjABAIxCMAEAjEIwAQCMQjABAIxCMAEAjEIwTVM4PKDS0iatWnVcpaVN6usbcLokAHAl\n25u4pqtAoFn19eskZailxZJ0VHV1DzpdFpAQ4b6wAvsDautvkz/Lr+qnq5Wdle10WUgTBNM0tbXN\nk5Rx/Sjj+jHgToH9AdUP1kszpZbBFmmfVPerOqfLQppgKm+a/P5hSWOPrrKuHwPu1NbfNv572Ogx\nkCSMmKapunq1pKNqa5snv3/4+jHgTv4s/+hIKUOSNXoMJAvBNE3Z2ZncU4rBQDis5kBA89raNOz3\na3V1tTKzuWdhquqnq6V9mnCPCUgWgglJ0RwIaF19/egX8JYWHZX0YF2dw1VhMtlZ2dxTgmO4x4Sk\nmNfWNv6Whea1cc8CwK0RTGkq2euyhv3+ca0jo8cAcCtM5blIODygQKB5QoNGdnbmLc9N9rqs1dXV\nOipNuMcEALdCMLmInbBJ9rqszOxs7ikBmBam8gwRj6k1O2HDuiwApmLEZIh4TK35/cPX/+3o4pPb\nhQ3rsgCYimAyRDym1uyEDeuyAHdx0/6GBJPDxhoWzp7tkvRDTWe0MxnCBkhfbtrfkGBy2LdTeIOS\nmlRQcFkPPzw35qk1Ox16AFJfW3+bNPP6QYrvb0gwOezbKbxMSQ9pyZLjqqtbE/Pr8pgOIL24aX9D\ngslhdhoW7OAxHcnlpvl9pCY37W9IMDksUd1xiQo83Jqb5veRmty0vyHB5LBENSzQDp5cbprfB5xG\nMLkUHXrJ5ab5fcBpBBNiRgegu+b3AacRTIgZHYDumt8HnMZeeYgZHYAA4inqEVN7e7sCgYAOHjyo\n2tpadXR0aHBwUOXl5RoZGdGePXuUk5OjxYsXa9OmTVGdk5ubG89rRYLQAQggnqIKpt7eXr399tua\nP3++RkZG9MEHH6impkbvv/++GhoadPnyZW3ZskXLli1TWVmZHnvsMVvnnDhxQg0NDSorK4v39SZN\nOt13oQMQQDxFFUz5+fl64YUX9NRTT6mvr0/5+fmSJJ/Pp1AopCtXrqioqEiSlJWVpYGBAVvnFBYW\nqqenJ+aLc5LJ913iHZp0AAKIp5ibHxYsWKBwOCxJ6urqks/n07Vr19TZ2Smfz6e+vj55vV7b53i9\n3kn/ZjAYjLXshIlEIgoGgzpzZmxqS5IydOaMZUTdkUhEpaXv6NChn2gsNMPht1RVVeJ0aTfeO1NR\nX2xMrs/k2iTz64u3mINpxowZWrlypXbt2qX+/n7t2rVLkUhEVVVVOnjwoNauXRv1OZMpLi6OteyE\nCQaDKi4uVklJu06f/va+S0lJRtR1x3OEEwwGdeFCtjT6oHNJw7pwIduI93TsvTMV9cXG5PpMrk0y\nv77W1ta4vl5MwbR//35J0ubNmyf83OPxaO/evRN+Fs05qSye913iPS146VKHpF9oLDQvXaqN+rUA\nIN5Yx5Qg8bzvEu927Ly8O/XFF9++Xl7enTG9HuyJx4avt3oNwC1Yx5QCRtuvretHsbdj3333yITX\nGz1Gsoxt+Noys0X1g/Xatm+bI68BmIoRUwqIdzs27d3OiseGr2waCzcjmFJAvNuxae92Vjw2fGXT\nWLgZwWSwdFqkm07iseHrrV7jQseFBFQLJB/BZDCTF+kievHY8PVWr0EwwS1ofnBQODyg0tImrVp1\nXKWlTerrG5jwezZHBZCOGDE5aKoREZujAkhHBJODphoR0T0HIB0RTA6aakRE9xyAdEQwOYgREQB8\nF8HkoGSNiMa3nXu9naqv/z+uazuPxzY/bq4HSCUEUxoY32QhWdq2zX1t52Nb9GimRhee7lPMLdlu\nqgdIJbSLp4F0aDtv628bf4mOb9FjWj1AKiGY0kC8N4E1kT/LP/4SHd+i53tzvzehnu/N+56D1QCp\nham8NDC+ycLr7VR19XqnS4q7eGzzE1dXJX0uabakK5J+4Gw5QCohmNLA+CaLYDDousYHKT7b/Ngx\nVXPD/1z5H+kuTTwGMC1M5QFRmOp5SKZNLQKphBETEIWpnodk3NQikEIIJiAKUz0PKdlTi4CbEExA\nFBgRAYlDMAFRYEQEJA7NDwAAoxBMaWggHFZTaamOr1qlptJSDfT1OV0SANzAVF4aag4EtK6+fvS+\nfUuLjkp6sK7O4aqAidgIN30RTEk0fpfvscdcOLHYdV5b2/ht3DSvjX3cTJTuH8xshJu+CKYkmupR\n6sky7PfLamkZ63TWsN8diz/d9kGe7h/MU60Vg3sRTElkyi7fq6urdVSjI6Vhv1+rq93R6uy2D/J0\n/2Ceaq0Y3ItgSqKpHqWeLJnZ2a68p+S2D/J0/2BmrVj6IpiSiEepJ5bbPsjT/YOZtWLpi2BKomQ9\nSj3RBsJhNQcCE6YCM7Odv5fjtg9yPpiRrggm2GZquzkf5IA7EEwxurkF/Nln850uKeFoN3c3t3U3\nIvUQTDG6uQU8HH5LjY0rnC4rodzabo5RbutuROohmGJ0cwt4e7v7v1kmo92cb+3OcVt3I1KPUcHU\n3d2tPXv2KCcnR4sXL9amTZucLmlKN7eAL1rk/n3nktFuzrf26Rkf4N5rXtVvr485wN3W3YjUY1Qw\nvfXWW9qyZYuWLVump59+Wj/96U81c+bMqf+hg25uAX/22SVOl+QKfGufnvEBrhnStn3bYg5wt3U3\nIvUYFUy9vb0qKiqSJGVlZWlgYEA5OTnfOS8YDCa7tNt68cVFN/47Eok4Xl9//5Cqqj5Xe3u2Fi3q\nU0XFEmVm3mFMfZMZX5v3mnd07/vr39q917yO123ie3em64w0tt1ixuhxPGp8cf2LN/77QscFXei4\nIEnqH+hX1f+vUnukXYvmLlLFoxXK9GRO+TvJzPdvjMm1SebXF29GBdPChQvV1dUln8+n/v5+ZWVl\n3fK84uLiJFc2fcFg0PH6SkubdOjQTyRl6PRpSzk5R1VXt8KY+iYzvrb67fXatm+bUfeYTHzvSgpL\ndHrw9I0ALyksSWiNpb8v1aFZh6RM6bR1WjnHc26M0G73O8nM92+MybVJ5tfX2toa19czKph+/OMf\nq6qqSh6PR2vXrtWMGan3uKj+/iGVljY5uoO4KXvyxYI1SdMzftrNe82b8Gm3202xMv2KeDEqmPLz\n87V3716ny4hJVdXnN0YrTu0gbsqefMmUrl184wM8GAwm/Jpv1xhB0wTixahgcoPRdvHkjFYme75T\nOu7JRxdfctyuMYKmCcQLwRRnixb16fTpxIxWbg6ikZEr+utfH9HNozO37MlnB9NIyXG7KVamXxEv\nBFOcVVQsUU5OYkYrN+8yUVBwSKl+LylemEYC3INgirPMzDtudMDF281NDZYV1uimQOlzL2kyTCMB\n7kEwpZCbmxrWrs3SnDnpdS9pMkwjAe5BMKWQ7zY1rE16KzoAJBrBlELSsakBQPpJvRWsAABXI5gA\nAEZhKg+Abem60waSg2ACYBs7bSCRmMpLAeHwgEpLm7Rq1XGVljapr2/A6ZJuaSAcVlNpqY6vWqWm\n0lIN9Ln/oYnpqq2/bfySOnbaQFwxYkoBN+/4cPnyQf3Xf813dAfzW2kOBLSuvn50lVVLi45KCX/S\nLZzBThtIJIIpBdy848OxY/3q7f1/cnIH81uZ19Y2/ku05rXxLdqt2GkDiUQwpYCbd3zIyMiRiXvk\nDfv9slpaxr5Ea9jPt2i3YqcNJBLBZNNkj5pIpJt3fBgZmam//tW8PfJWV1frqEZHSsN+v1ZX8y0a\ngH0Ek0033+9JxjTazTs+9PUNGLlHXmZ2NveUAMSMYLLJhMeWszXRtwbCYTUHAhNGaZnZrKcBUhnB\nZFM6PrbcZHQCAu5DMNmUSo8tv9Vowm3oBATch2CyKZWm0W41mlj04otOl2XLVFN1sXQCmrCtjgk1\nAKYhmFzMDaOJqabqYukENGFbHRNqgDul8pcegsnF3LCuaKpwjaUTsK2/TZr57Ys7sa2OCTVEK5U/\n+NJBKn/pIZhc7FajifMXLsT0msnugktkuJqwrY4JNUQrlT/40kEqf+khmFzslqOJGIMp2V1wiVy0\na8K2OibUEK1U/uBLB6n8pYdggi3Jvm813am6aEZyJmyrY0IN0UrlD750kMpfeggm2GLqfSvWMyVf\nKn/wpYNU/tJDMMEWU/fDc0MHYqpJ5Q8+mI1ggi2m7odn6kgOgH0Ek+Gc2M08FZk6kgNgH8FkOCd2\nM09Fpo7kANg3w+kCcHsm7GYOAMlEMBludPdy6/oRu5kDcD+m8gyXSruZA0A8EEyGS6XdzAEgHqIK\npsbGRn388ceqrKzU119/rZ07dyo3N1eZmZl6/vnndejQIZ08eVKRSETPPfecCgoKbJ3j8Xi0devW\neF8rJkHnHwCT2A6mI0eOqKenR1evXr1x/MADD2j9+vWqqKhQd3e3GhsbdeDAAZ0/f1779u1TSUmJ\nrXPKy8sVCoXk9XrjfsH4Ljr/AJhkymBqaGjQ4cOHJUkZGRmqqKjQvffeq/3790uSent7tWLFCklS\nQUGBenp6NHv2bElSYWGhQqGQLl68qOXLl0/7HK/Xq56enkmDKRgMxnLNCRWJRFKuvjNnxh4VL0kZ\nOnPGcuQaUvG9Mwn1Rc/k2iTz64u3KYNp48aN2rhx44SfdXR03PjvoqIidXd3S5JCoZB8Pp9mzBht\n9uvs7JTX61VRUZG6urpsnXO70VJxcbGda0yqYDCYcvWVlLTr9OmxcLJUUpLhyDWk4ntnEuqLnsm1\nSebX19raGtfXi7n5Yd26daqsrNSHH36oRYsWqaCgQI8++qheeuklDQ0Nafv27crLy1NlZaVOnTpl\n6xwkB51/AEySYVmWNfVp5mhtbb0x5Wci07/ZmFyfybVJ1Bcrk+szuTbJ/Pri/bnMAlsAgFEIJgCA\nUQgmAIBRCCYAgFEIJgCAUQgmAIBRCCYAgFEIJgCAUQgmTMtAOKym0lIdX7VKTaWlGujrc7okAC7F\n85gwLc2BgNbV14/uptfSoqOSHqyrc7gqAG5EMGFSA+GwmgMBzWtrU+Tzz8ftPy7Na2tzsjQALkYw\nYVLjR0lHJH27/7g07Pc7WhsA9yKYMKl5bW03Rkn/V1J9fr6+f9ddGvb7tbq62snSALgYwYQJxk/f\nnQ2F9EONjpI8kgrXr9ca7isBSDCCyUXC4QEFAs0TnquUnZ1p6zXGT9/9t6TaxYt1l9fLKAlA0hBM\nLhIINKu+fp2kDLW0WJKOqq7uQVuvMX76LkvSXV6v1rz3XpwrBYDJsY7JRdra5knjeudGj+0Z9vs1\n9uRImhwAOIERk0uEwwMKhc5KN+4KWfL7h22/zurqah3V6MiJ6TsATiCYXCIQaNYXX2yUdFTSXC1e\n/Lmqq39i+3Uys7ONWTg7vhFjLCQzs7OdLgtAghFMLjE6bZclafSektebYbvxwTTsNgGkJ+4xucTo\ntN23d4eimcYzzfhGjFTabYJ9BYHYMGIySCzt3tXVqyUdnfBvU92w3y+rpSXldptgpAfEhmAySCzt\n3tnZmbZbw02Xqo0YqTrSA0xBMBkkHu3ebmJSI4YdqTrSA0xBMBnE7x++PlKKvt0bzkvVkR5gCoLJ\nIG68T5SOUnWkB5iCYDKIG+8TAYBdtIsDAIxCMAEAjEIwAQCMQjABAIxCMCFtsFUQkBroykPaYKsg\nIDUwYkLaYKsgIDUQTEgbPJ0XSA22p/JeeeUVhcNhffXVV3ryySe1ZMkS7dy5U7m5ucrMzNTzzz+v\nQ4cO6eTJk4pEInruuedUUFBg6xyPx6OtW7cm4nqRxtgqCEgNtoPpzjvv1MMPP6xPP/1UR48e1dmz\nZ/XAAw9o/fr1qqioUHd3txobG3XgwAGdP39e+/btU0lJia1zysvLFQqF5PV6E3HNSFNsFQSkhimD\nqaGhQYcPH5YkZWRkqKKiQhcvXlRtba127NihgwcPasWKFZKkgoIC9fT0aPbs2ZKkwsJChUIhXbx4\nUcuXL5/2OV6vVz09PZMGUzAYjPGyEycSiVBflEyuTaK+WJlcn8m1SebXF29TBtPGjRu1cePGG8et\nra2qr6/Xyy+/rJycHBUVFam7u1uSFAqF5PP5NGPG6K2rzs5Oeb1eFRUVqaury9Y5txstFRcXR3m5\niRcMBqkvSibXJlFfrEyuz+TaJPPra21tjevr2ZrKGxwc1NatW3X//fdr9+7dWrlypR555BFVVlbq\nww8/1KJFi1RQUKBHH31UL730koaGhrR9+3bl5eWpsrJSp06dsnUOJhoIh9UcCEy4R5KZne10WQAQ\nV7aCyePxqLm5+Ts//93vfjfh+KGHHtJDDz0U8zmYKBnrcAg/AE5jgW0KScY6HBahAnAa65hSSDLW\n4bAIFYDTGDGlkGSswxn2+2W1tFx/uDuLUAEkH8GUQpKxDodFqACcRjBhAhahAnAa95gAAEYhmAAA\nRiGYAABGIZgAAEYhmAAARiGYAABGIZgAAEYhmAAARiGYAABGIZgAAEYhmAAARiGYAABGIZgAAEYh\nmAAARiGYAABGIZgAAEYhmAAARiGYAABGIZgAAEYhmAAARiGYAABGIZgAAEYhmAAARiGYAABGIZgA\nAEYhmAAARiGYAABGIZgAAEYhmAAARiGYAABGmWX3H/zlL3/R2bNnNTw8rG3btik7O1s7d+5Ubm6u\nMjMz9fzzz+vQoUM6efKkIpGInnvuORUUFNg6x+PxaOvWrYm4XgCA4WwH04IFC/TEE0+oqalJLS0t\nsixLDzzwgNavX6+Kigp1d3ersbFRBw4c0Pnz57Vv3z6VlJTYOqe8vFyhUEherzcR1wwAMNiUwdTQ\n0KDDhw9LkjIyMlRRUaF//vOfevXVV/XHP/5Rx44d04oVKyRJBQUF6unp0ezZsyVJhYWFCoVCunjx\nopYvXz7tc7xer3p6eiYNptbW1hgvO7GoL3om1yZRX6xMrs/k2iTz64unKYNp48aN2rhx443j9957\nT2vXrtX999+vHTt2aN26derq6pIkhUIh+Xw+zZgxeuuqs7NTXq9XRUVFts+ZLJTGwgsA4E4ZlmVZ\ndv7BK6+8olAopG+++Ubr16/XihUrVFlZqaysLOXn5+uZZ57RkSNH1NzcrKGhIW3fvl15eXlRnQMA\nSD+2gwkAgESiXRwAYBSCCQBgFNvt4olg+tqoV155ReFwWF999ZWefPJJLVmyxKj6Ghsb9fHHH6uy\nslJff/21UbVNV3d3t/bs2aOcnBwtXrxYmzZtSujfk6T29nYFAgEdPHhQtbW16ujo0ODgoMrLyzUy\nMvKdeqI5Jzc3N6raTp06pTfffFMej0d5eXmaO3euMfW1t7erurpaeXl5Wrp0qS5dumRMbeP96le/\n0o9+9CN1dnYaVd+FCxf07LPPqri4WF6vV1lZWcbU19HRoVdffVUej0dZWVnO/e/OMkBTU5NlWZZ1\n5MgR629/+5v19ttvW3//+98ty7Ks8vJyq6uryyotLbUsy7L+85//WJWVlbbP2bFjh9Xd3R1Vfe+8\n845lWZb1ySefWHv37jWqvn/84x/Wvn37rJ07d1qWZRlVmx3V1dXWqVOnLMuyrKeeesr65ptvEvr3\nenp6rN///vfWE088YV2+fNkqKyuzLMuyTpw4YdXU1Eyo5+mnn7YGBwdtnfPvf//bqqmpibq+f/3r\nX9bQ0JBlWZb185//3HrmmWeMqe+TTz6xzp8/b125csW42sYcOHDAevnll6133nnHuPoaGxutX/7y\nl9aOHTusd99916j6fvOb31g1NTXWjh07rOPHjztWmyNTeQ0NDdq8ebM2b96sLVu2aNGiRTfWRpWU\nlKi3t1eFhYWSbr/uyc45Y2ujoqnvzjvv1MWLF1VbW6uf/exnjtZ3c23f//739fDDD9/4vdPvXbR6\ne3tVVFQkScrKytLAwEBC/15+fr5eeOEFzZ8/X319fcrPz5ck+Xw+hUKhW9Zj55zCwsKY3rM1a9Zo\n/vz5qqmp0YYNG7RgwQJj6lu6dKlmzZqlsrIyrVy50qjaJOndd99VZmamli1bpmvXrhlX37333quq\nqipVVVWprq7OqPra29u1Zs0a/fa3v9Xrr7/uWG2OTOWZtjZqqvpaW1tVX1+vl19+WTk5OSoqKlJ3\nd7cj9d1cmzQ6/B7jZG2xWLhwobq6uuTz+dTf36+srKyE/r3xFixYoHA4LEk3arh27Zo6Ozvl8/nU\n19cnr9dr+5xY3rOhoSHt3r1bGzZs0H333adjx44ZU99nn32mhQsXqra2Vlu3btXVq1eNqU2SDh8+\nrOzsbJ07d06SdMcddxhVXzAY1LJlyyRJc+fOVSgUMqa+goIC3XHHHZo5c6bmzp3r2P8vjGgXN3lt\n1ODgoNatW6f7779fM2bM0MqVK/XII48YU580Gkz79+/Xr3/9a12+fNmo2qart7dXVVVV8ng8uuee\ne/T4448n9O+Neeqpp7R//3698cYbOnfunPr7+7Vr1y5FIpHv1BPNOR6PJ6q6Kioq9OWXX2rhwoWa\nOXOm7r77bmPq++ijj/T666+rsLBQs2fPVmFhoTG1jdfY2Kg5c+bo4sWLRtX36aefav/+/crPz9cP\nfvADRSIRY+r74osv9Oc//1mZmZlauXKlvvrqK0dqMyKYAAAYQ7s4AMAoBBMAwCgEEwDAKAQTAMAo\nBBMAwCgEEwDAKAQTAMAo/wv2oI+pffJ28AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbc8c8b6f28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(cbc5_params[cbc5_prediction==0,0],cbc5_params[cbc5_prediction==0,2],c='b')\n",
    "plt.scatter(cbc5_params[cbc5_prediction==1,0],cbc5_params[cbc5_prediction==1,2],c='g')\n",
    "plt.scatter(cbc5_params[cbc5_prediction==2,0],cbc5_params[cbc5_prediction==2,2],c='r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def eval_IPL_overlap(prediction):\n",
    "    current_overlap_IPL=np.zeros(69)\n",
    "    for i in range(69):\n",
    "        current_overlap_IPL[i]=np.sum(overlaps_IPL[((overlaps_IPL[:,0]==i+390+144)&(prediction[overlaps_IPL[:,1].astype(int)-390-144]==prediction[i]))\\\n",
    "                        |((overlaps_IPL[:,1]==i+390+144)&(prediction[overlaps_IPL[:,0].astype(int)-390-144]==prediction[i])),2])\n",
    "    current_overlap_IPL/=cbc5_ipl_area\n",
    "    return(current_overlap_IPL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def eval_OPL_overlap(prediction):\n",
    "    current_overlap_OPL=np.zeros(69)\n",
    "    for i in range(69):\n",
    "        current_overlap_OPL[i]=np.sum(overlaps_OPL[((overlaps_OPL[:,0]==i+390+144)&(prediction[overlaps_OPL[:,1].astype(int)-390-144]==prediction[i]))\\\n",
    "                        |((overlaps_OPL[:,1]==i+390+144)&(prediction[overlaps_OPL[:,0].astype(int)-390-144]==prediction[i])),2])\n",
    "    current_overlap_OPL/=cbc5_opl_area\n",
    "    return(current_overlap_OPL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def eval_clustering_GMM(prediction):\n",
    "    distances_to_means=np.zeros(69)\n",
    "    for i in range(69):\n",
    "        distances_to_means[i]=np.sqrt(np.dot(cbc5_params[i,:]-cbc5_GMM.means_[prediction[i],:],\\\n",
    "                                     np.dot(np.diag(1/cbc5_GMM.covars_[prediction[i],:]),cbc5_params[i,:]-cbc5_GMM.means_[prediction[i],:])))\n",
    "    return(np.sum(distances_to_means))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "l1=1;l2=1.5;l3=1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cbc5_ipl_area,overlaps_IPL=calculate_overlap((BC_IDs[:,3]>62)&(BC_IDs[:,3]<66),'IPL')\n",
    "cbc5_opl_area,overlaps_OPL=calculate_overlap((BC_IDs[:,3]>62)&(BC_IDs[:,3]<66),'OPL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cbc5_overlap_opl=np.concatenate((overlaps_OPL,overlaps_OPL[:,[1,0,2]]))\n",
    "cbc5_overlap_ipl=np.concatenate((overlaps_IPL,overlaps_IPL[:,[1,0,2]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iterations: 174\n",
      "switches: 13\n"
     ]
    }
   ],
   "source": [
    "num_iter=0\n",
    "num_switches=0\n",
    "switch_old=[]\n",
    "\n",
    "while True:\n",
    "    current_IPL_overlap=eval_IPL_overlap(cbc5_prediction)\n",
    "    current_OPL_overlap=eval_OPL_overlap(cbc5_prediction)\n",
    "    cost=l1*eval_clustering_GMM(cbc5_prediction)+l2*np.sum(current_IPL_overlap)+l3*np.sum(current_OPL_overlap)\n",
    "    switch_array=np.argsort(current_IPL_overlap+current_OPL_overlap)\n",
    "    switch_array=switch_array[np.in1d(switch_array,switch_old,invert=True)]\n",
    "    switch=switch_array[-1]\n",
    "    prediction_new=np.vstack((cbc5_prediction,cbc5_prediction))\n",
    "    sorting=np.argsort(cbc5_GMM.score_samples(cbc5_params[switch,:].reshape(1,-1))[1][0])\n",
    "    prediction_new[:,switch]=sorting[sorting!=cbc5_prediction[switch]]\n",
    "    if l1*eval_clustering_GMM(prediction_new[0,:])+l2*np.sum(eval_IPL_overlap(prediction_new[0,:]))+l3*np.sum(eval_OPL_overlap(prediction_new[0,:]))<cost:\n",
    "        cbc5_prediction=np.copy(prediction_new[0,:])\n",
    "        switch_old=[]\n",
    "        num_switches+=1\n",
    "    elif l1*eval_clustering_GMM(prediction_new[1,:])+l2*np.sum(eval_IPL_overlap(prediction_new[1,:]))+l3*np.sum(eval_OPL_overlap(prediction_new[1,:]))<cost:\n",
    "        cbc5_prediction=np.copy(prediction_new[1,:])\n",
    "        switch_old=[]\n",
    "        num_switches+=1\n",
    "    else:\n",
    "        new_overlaps=cbc5_overlap_ipl[(cbc5_overlap_ipl[:,0]==switch+390+144)&np.in1d(cbc5_overlap_ipl[:,1],np.where(cbc5_prediction==prediction_new[0,switch])[0]+390+144)]\n",
    "        new_overlaps[:,2]+=cbc5_overlap_opl[(cbc5_overlap_opl[:,0]==switch+390+144)&np.in1d(cbc5_overlap_opl[:,1],np.where(cbc5_prediction==prediction_new[0,switch])[0]+390+144),2]\n",
    "        new_overlaps[new_overlaps[:,2]==np.max(new_overlaps[:,2]),1][0]\n",
    "        prediction_swap=prediction_new[0,:]\n",
    "        prediction_swap[new_overlaps[new_overlaps[:,2]==np.max(new_overlaps[:,2]),1][0].astype(int)-390-144]=cbc5_prediction[switch]\n",
    "        if l1*eval_clustering_GMM(prediction_swap)+l2*np.sum(eval_IPL_overlap(prediction_swap))+l3*np.sum(eval_OPL_overlap(prediction_swap))<cost:\n",
    "            cbc5_prediction=np.copy(prediction_swap)\n",
    "            switch_old=[]\n",
    "            num_switches+=1\n",
    "        else:\n",
    "            new_overlaps=cbc5_overlap_ipl[(cbc5_overlap_ipl[:,0]==switch+390+144)&np.in1d(cbc5_overlap_ipl[:,1],np.where(cbc5_prediction==prediction_new[1,switch])[0]+390+144)]\n",
    "            new_overlaps[:,2]+=cbc5_overlap_opl[(cbc5_overlap_opl[:,0]==switch+390+144)&np.in1d(cbc5_overlap_opl[:,1],np.where(cbc5_prediction==prediction_new[1,switch])[0]+390+144),2]\n",
    "            new_overlaps[new_overlaps[:,2]==np.max(new_overlaps[:,2]),1][0]\n",
    "            prediction_swap=prediction_new[1,:]\n",
    "            prediction_swap[new_overlaps[new_overlaps[:,2]==np.max(new_overlaps[:,2]),1][0].astype(int)-390-144]=cbc5_prediction[switch]\n",
    "            if l1*eval_clustering_GMM(prediction_swap)+l2*np.sum(eval_IPL_overlap(prediction_swap))+l3*np.sum(eval_OPL_overlap(prediction_swap))<cost:\n",
    "                cbc5_prediction=np.copy(prediction_swap)\n",
    "                switch_old=[]\n",
    "                num_switches+=1\n",
    "            else:\n",
    "                switch_old.append(switch)\n",
    "    num_iter+=1\n",
    "    if len(switch_old)==cbc5_prediction.shape[0]:\n",
    "        break\n",
    "print('iterations:',num_iter)\n",
    "print('switches:',num_switches)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "CBC5T=BC_IDs[(BC_IDs[:,3]>62)&(BC_IDs[:,3]<66),0][cbc5_prediction==1]\n",
    "CBC5O=BC_IDs[(BC_IDs[:,3]>62)&(BC_IDs[:,3]<66),0][cbc5_prediction==2]\n",
    "CBC5I=BC_IDs[(BC_IDs[:,3]>62)&(BC_IDs[:,3]<66),0][cbc5_prediction==0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "BC_IDs[np.in1d(BC_IDs[:,0],CBC5T),4]=63\n",
    "BC_IDs[np.in1d(BC_IDs[:,0],CBC5O),4]=64\n",
    "BC_IDs[np.in1d(BC_IDs[:,0],CBC5I),4]=65"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# np.savetxt('data/BC_IDs_new',BC_IDs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot density distributions of the axon terminals of the new CBC5 subtypes (Supp. figure 7B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN0AAACACAYAAACP1GwzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXd8VFXexr93espMkklvpEJIQgcpgiBFmmIDxQor6qpr\nYRUsILo2ZIFXXcWCK4riLrJ2RVFEBBRCDSVAEgIJBEjvk5lk6r3vH4GslEDKTBLc+/18YJKZe899\n7uQ+95x7zu/8jiBJkoSMjEy7oehoATIy/2vIppORaWdk08nItDOy6WRk2hnZdDIy7YxsOhmZdkbV\n0h3S09M9oUNG5pKmf//+zd62xaZr6QHak6ysLJKTkztaRpN0Zn2dWRt0bn0trYjk5qWMTDsjm05G\npp1pVfNSpvm4XC5eeeUVTCYTpaWlJCUlMWvWLObMmUNdXR0ajYaqqipefvllQkJCePfddzl27BiS\nJKHRaHj22WdRqRr+THa7ncmTJ9OjRw/UajUzZ87k73//OxUVFVRVVdG1a1emTp3aaZv/Mg3IpjtF\nUbkFS72j1fv7eKkJD/I55/1Vq1aRkpLCNddcQ1ZWFkeOHKGqqgqAZ599lsDAQD744AN27dqFt7c3\nAAsWLABg06ZNFBYW0qVLFwAyMjLw9fVFqVTSu3dvAgMDWbx4MTt37mTv3r3ce++9rdYv037IpgNq\nzDbu//vPiG0I/VYIsOK58fj5as94Pycnh3vuuafx90mTJjX+PH/+fJRKJfv37+f999/n+++/p1ev\nXo2fjxgx4oyyoqKiWLRoEdHR0Tz66KMMGzaM0NBQ5Jj1SwvZdICfr5alT41pc013tuEAkpOT2bJl\nC7fccgsAH374IcOHDwfg6aefJjAwkPT0dN58800mTpxIWloaQ4YMAWD16tWEh4czYMAAAPLy8vDz\n8wPAYDDgcrlarVem45BNd4rzNQ3dwZQpU5g/fz6zZs3CZDKRkpJCfHw8AC+88AI6nY6amhqmT5/O\nkCFD2LVrF48++ihqtRovLy+eeeaZxrJiYmJYuHAh4eHhhIeHExER4RHNMp5FaOl8uvT09E77oN6Z\nx3Kgc+vrzNqgc+trqSfkIQMZmXZGNp2MTDsjm05Gpp2RTSdzyWOpOUFt5ZGOltFs5N5LmUsSp6MO\nh7WGypJ9FOetByAq6VpCY67oYGUXRzadh/FkGNgLL7xAYWEhr7/+OhqNhpqaGgYNGsTtt9/ewWft\nWSw1Jzic/k9cTisAAWG9cdjMnDz0LTrvIPyCO2cv52lk052i2FxGnb2u1ft7a7wJ8w0+531PhoGd\n3vb555/HaDQC8Pzzz3fqYZ224nTUc2T3+2i8AgnpMhR9QDxa70AkSeJw+j/Jz/qSFL+/otJ4ZtzV\nHcimA0zWWmau+VubwqkEQeC9axdi0OnPeN+TYWBDhw6lrq6u0XAAI0eOJDMz8w9ruvKC7bicVpL7\n3oVG59f4viAIxKTeRNa21zm6fyWJ/e5BEIQOVNo0sukAg07P6xOfb3NNd7bhwHNhYH5+frhcLjQa\nDRUVFQQGBgKwefNmJkyY0Orz6OxUFWdgCEo6w3Cn0XoZiU29mdy9H1JbeRhDYLcOUHhxZNOd4nxN\nQ3fgqTCwsLAwIiMjmTt3LvPnz0en02Gz2ejXrx99+/b1yLl0NLb6SupMJ4jtcUuT2/gFp6DzCaXs\nxLZOazo5DKwd6cz6OrM2aNDnIxyisnA3PYfPRanSNbltUe46io9tpM/IFxAUSo9rk8PAZP6QSJJI\nVfE+AiP6X9BwAL4BCYguO/XmonZS1zJk08lcEoi2Ypx2MwGhvS66rY9fFAgKzNX57aCs5cimk7kk\ncFpy0HoH4eMfe9FtFUoN3voILLLpZGRahyS6cNWfwBjet9nDAL7+sZhrOqfp5N5LD+PpiJQ5c+Yw\ne/bsxiGDPyLm6nyQHPgFdm/2Pj7+MZQe34zDZkKtNXhQXcuRTXeK+qJiXBZLq/dX+vjgFR52zvue\njkj5X6C6dD8oNHj7RTV7H1//OABqq/IwhvXxlLRWIZsOcNTUsPsvD4Motr4QhYKBHy5D7XfmoK2n\nI1L+6DhsZspObEXlm4IgNP9pSKPzQ+sdRG1lrmy6zojaz49+by9pc013tuHAsxEpYltuEpcIZSe3\ngiCg1vds8b76gARqK3M9oKptyKY7xfmahu7AkxEpf/TERJLoouzEVgLD+1OnuPDY3PnQGxMoL9iO\n3Vpz3rCxDkNqIbt27WrpLu1GZmZmR0u4IJ1ZX2fUZqo4Iu1aO1syV+W3Sp/dWiPtWjtbqijc7QF1\n/6WlnrgkhgwkUWJX2jHef2Mz33+egd3mpLighoxdJxBdf/wm1v8qpooclGrvFnWg/B611oDWO7jT\nzSq/JJqXv60/zMYfDxESpicj/SSZ+wqx21y4XCKFJ2sYf32PjpYo4wFMFTkYjF1b1IFyNnpjArUV\nh92oqu10etOVFteyce0hBo+IZ+y1qRw9XM5P3xwkMTkEpUrBb+ty6N4zjNiEoI6WKuNGnHYLdaYC\ngqOGtKkcvTGR8pPbsFur0ej83aSubXRq06VvzSdtwxECjN6MmtAwMBrXNYj7Zjd0pTudLvJzK/jP\nBzt58KlRHSlVxs2YKo8AEobArm0qRx/Q0GlVW5lLYETnmB3TaZ/pDh0o5vsvMvA3ejNl2gBU6nOn\naKhUSm6a1h9BEPhlTVYHqJQ5G1NWNvvnzCPn1dexlpa2vpzyQ2i9g9F4BbRJj1qrR+cT2qmGDjpd\nTedyinz60S4OZ5YQmxjEbfcMQqlq+t7g7avlynFJ/PjNAQLDo6ATTglbsmQJhYWFlJeXExYWxhNP\nPIHJZGLGjBmNkSUpKSmkpKSwYMECunZtuLsPHTqU6667rrGcBQsWYLFYcLlcpKamcscddwDwySef\n8OWXX/LZZ58BZ4aeWa1WwsPDmTVrlsfPs2b/AQ698g9cFgvmvKNU791H0pOz8UtNaVE5kiRSU55F\nYLh7aia9MYGaskwkSeoUKRw6nen27y7gcGYJ19/ah579o5r1JQ24PIb0bfnsSStj8FARpbLlFXhl\nuQVrG1bt0XmpMZ5nEZKvv/6asLAwHn74YbKysrDZbCxcuJAHHniAIUOG8NxzzzVuu2PHDiZOnHje\ndebq6+s5evQob731Fmq1mi+//LLxszVr1nDVVVexdu1axo0bd0boGTQMsldXV+Pv77lnmvLNW8h5\n7Q18uyaS9PgsBKWSrBfnk71gEf2XvonK17fZZVnNJTjt5jY3LU/jF5xM2Yk06msL8TZEuqXMttCp\nTOdyimzblEt8t2B6DYhu9n4KpYKJN/bko7fTePeVTdz8p8sICmn+H9litvHW33+hLcu8CQI89txY\nfM5aLisjI4O77rqr8fc+ffrw5ptvArBt2zbmzJkDwIwZMwD44YcfyMvLA2DWrFkEBTV0EHl5efHI\nI4/w2muvYbPZGgfYN23aRO/evbnxxht56KGHGDdu3AVDzzxB3fETHHrlHwQOGUy3Rx9BoVYD0H3u\nk+y+/yFOfv4lsX+a1uzyaitzEQQlPqfiJ9uK3piIUqWjqiRDNt3Z/PbzYcpKzUya2vJYuZiEQEZd\nG8WeLZV8+XE6d8+84oLN0t/j46vlwadGtbmmO9tw0NBs3Lp1K9HRDTeRrKysxkiSwYMHn1PTTZgw\n4bw1XXZ2Njk5OTzxxBMAPP7444wePZoVK1bg7+/PkiVLqKysZOvWrecNPRsxYgRxce65iM/m+MpV\naIzGMwwHoA0MJOL6ayn44iuChg/D99SN4mLUVh3Bx68LSpXGLfoUChV+walUlewnInF8hzcxO43p\n7DYn237NY9AVcUR2aV0zKDBEx4139Gf5ks1s2XCE4Vc1PzHN+ZqG7mDy5Mm89tprPP7445jNZoxG\nI3PmzKGmpqZF5SQmJvLpp5+yefNmVCoVer0ek8mEj48Pr7zyCgCHDh1i8eLFvPvuu7z00kvMmjUL\nQRCIjIz0mOFMmVlUbN1GwoP3n2G400TfPIWKLWnkr/g3qc89c54SzkSSRGor8wjp4t5g7oDQXlQW\npWM1l+Cl90zIX7PxdMhLczm4t0B6/rFvpYoyc6vLOB0qtOaLDGnh0z9IdpvDXfLcQmcMtTpNa7Ud\nevUf0s5775dEl6vJbcp+2yxtvvZGqTxt60XLqzeXSrvWzpaqy7Lcou80Lqdd2v3z01LB4bVtKud8\nXLJhYIezSgkO9XVLjTNoeDzWegcH9hS6QZlMUzhqaqhI20bomNEIiqYvpcChlxMwoD/Zf1/MgXl/\nw2WzNbmtxXQCAB9D85/pm4NCqcYvqDs15R0/tNQpTOdyiRzJKiGhe4hbyjMG+ZDQPZgdm4+2KWuz\nzIUpT9uK6HQSOvaqC24nCAJd//owXW6/ldpDORxd9kGT21pqjqPRBXgkLbpvQDx1tYW4nHa3l90S\nOoXpDh0oxmK207Of+3qWBg6Lo6TQxPG8SreVKXMmFVu24tcjFY3/xafNqPV6om+eQvStUyldvwF7\ndfV5t6utOIzemOhuqQD4BsSCJFJnOu6R8ptLpzDdzi3HiIoJIDzKfeNIiUkhGIN82LH5qNvKlPkv\ntopKag5mEjT08hbtFzpmNCgUFP/40zmfOexmrJZS9MYEd8k8Ay/fMBQqHeaqjr0mOrz3sqTIRH5u\nBTfc5t5U4IJCoP/lMfz8XRY2qxOtruNOtamIlMmTJ/PNN98QGhrKmjVrsNvtXH/99Y37deYIlJJ1\nP6NQqQga1jLTqQ16Qq8aTcFX3xA4eBA+sTGNn9WbCgDw8XPv89xpBEGBr18XzNXHPFJ+c+lQ0zns\nTjb+kI3eoCOlt/tnQSd0C2admMmJY5UkXuR50VZXjvPUemetQaXSofU+d6bDhSJSUlJSmDdvHu+8\n8845+3XmCBTR6aRk7TqCRwxvUaTJaWLuvB1TZhYH5j1L9zlPNoaJ1dUWoFBqzvs9ugvfgDiKj25E\nksQ2TRlqCx1mulqTlX++sgmL2c6km3s3eyC7JQSH6vHyVpOfV3FB0znsZg5sXgS0pdNFoNeVz6LW\nnHkRNhWRIggC0dHRXHHFFSxatIg+fc4MCOjMESiV23dir6wkbOK4Vu2v8vamx4vPkzV/AYcWv0q/\nt5eg8vaizlSIlz7Co2bw9Y9FdNmoNxfjre+YdBcdZrqfvjmIyyVx10NDiY4zXnyHViAoBGISAjme\nW3HB7dQaX3oMe6LNNd3ZhoOmI1KkU8G3Y8aMYffu3fz444+MGvXf6UmdIQLF5rSzJucXNh/fSf+I\nngyM7INDdOD84Qf0SUlnRJiYbGbyq0+SGtwNxQWGD06jNujp9thM9jz0V05+9jmx0++krrYAQ1BS\nmzRfDG9Dl4aU61VH/7dMV11Zx8F9hUy4oafHDHeamPhA1n2XicPuRK1p+nS13kGcG8TVdpoTkTJr\n1iz+9Kc/nbFfR0egmGxmXtz4OserC+gWGMfq7HV8nbUW73oX9+yv4MDoeD784QVCfIOIMoSz9vBG\nbC47fjoDAyN7E+gdgFqhZmTcEHy15+/+14WGEnnDdZz84isCrxyKra4cb71n50UqVRq89ZGYq4+5\nPeqluXTIUlkb1x5i68ZcHvvbWLd2cJxvuafighr++eqv3Hn/EOK6duzs8s68HNVpbaIocqTyGN9m\nr2N/aTbPj3yM2IBoqupryKnIw/rbTpSr1pLx8Bg0/v5klR3mhKmIXqHJXBEzkIySLPaXZFNtNQGg\nUqhINMZwVcJw+kf2xFvtdcZxXTYbux94CN8h3amLKSZ58F/PG5Tszu/u5KHVVBSm02vEM25ZSqul\nnmj3mk50iezdfpzU3hHt0qMYEm5A56UmP7eiw013KfDP9JX8krcFjVLNPf1vJTagoVkc4OXHoKi+\nZOWtxZGUxF/GPACAJEnUWE346QwIgsCIuMEAWB1WCmtL2FmQwReZa8guz0WlUDGx20im9piEWtkQ\np6nUagkdM5qiI7+gjg3Ey9fzcZHG8H6U5P9KdVkmAaEtz6fZVtrddIezSzHVWOk3JObiG7sBhUKg\nS5yRY7nlgGefFzor9TYneQU1qJQCJZV1nCgxExbozagB0WdE3O8uPMAveVu4tvtV3JA8Hh+N9xnl\nuGw2qvfuI3rqTY3vCYKAv9e5g+M6tY54Ywzxxhim9pzEsaqTbDiaxrfZ6yiqLWXW0D+jONVhEjLq\nSoorNqLGt10WcfQ2ROJtiKL85Pb/DdPt3ppPaISh1TMJWkN8UjA/fXMQa70Dnde5kfB/ZPKLTbyw\nbBulVfWN73lpldTbXOzMKuGRm/vgrVPjkkTe372KnqFJ3N7rhvNOf6nJ2I9ot2MceFmLdcQGRHFX\nwM0kByfyatp7/JCzgauTRgOgDQ1FGeWLVOFs/Ym2kKDIgRzP+gp7fVWbU0K0lHYdqKipquNIdin9\nh8S065ymrsmhiKJE7qGydjtmZ6DSZOXZd9PQaVUsfGgYr8wczr+eH8+nL1/Dk9MGsONgMQ8u3kDO\n8SqO1OZTZqng1p7XN/m3qdyxE11YKF7RrctDCTA4uh9j4ofx6cHvOFHTEJBeby4CHwHrvgKkdkoV\nbwzvi0KhorxgZ7sc7/e0q+n2bD+BSq10a4xlcwgI9CY41JfDmSXtetyOZtW6QzicEi/8eQgpcYF0\n6xKA36mJtsN6R/LW46MINOiYt3QLm4syiPGPIsF4/ma/JIpU7txFwGWXtfmGeVvv6wnyNvLSxjco\nt1RSVZKBQlDjyK7AnJvXprKbi1KlIyC8D+WFO5Gk9k1Y3G6mkySJ/btPktI7Aq2u/Zt43XqEcehg\nMXZb+zVhOhJznZ0Nu04wfkgMgX5e590mPMiHF++/nPBwJcet+QwMa9pQ5tw8HFXVGAcOaLM2X40P\n80Y8glKh5MWNr1NWuAe/kFSUWi9yf9nCpz/nsHV/ocdniARFDsRhrcZUfsijxzmbdjNd4Ylqqirq\n6NG3YwYk+w+OwW5zsm/niQ45fnvz47Z8nC6Rq4deeLzOS6ti4DAnSAoKDjU9W6Byx06UPt4YUtzT\nbR/g5ce8EQ8TIIDLWsmRGiPZ2nCOr9vIJ2uzefnDnTz0fxtYuy2f0qq6RgPaHC7yCmpwuiGdvo9f\nDF6+4ZTk/9bmslpCu3WkHNhdgI+vhrjEjum29zd6071nONt/O8qAy2MRFB2fis1TuESJNWlHuaJP\nZJO13GmKa0v5Ke8XwoV4Nu4s4dYxdYQavc/ZrnLHTgL690Ohcs8lY7U52bLTRHh5CPaQoyw7uJ2B\niX3ouutLlt0Rz1FNMN9tzuPNz/YC4K1V0DOxlowjZdTbXBh8NNw8phvXXhHf6uauIAiExY3k6P6V\n1NcWt1sah3ap6USXyIE9BaT2jUTRivR47mLQFXFUllvI2H2ywzS0BxmHyyirqmfc4NgLbidJEm9u\n/wiD1pdbk0bi46Xms/U552xnLS2l7lg+xsta3mt5Pqw2J/OWprFybTaJAdXYtGFgLCHhhkC0oSGU\n/7yOAcmhPHfvEN6bO4ZnZgyif1c/TBY744fE8ezdgxiUGsaybw6w4KOdWO2tf2TwD+2JUqWjsniP\nW86tObSLA/IOl5+apNr6Xi930CU+kJTeEfz0zUHqzE2nDLjU+XnncSKCfEi5SIhdRkkWORV5zOg/\nFT+dFzdemcjPO45TUll3xnaVO3YhKJUE9Gv79CuHU+S5Zds4Vmzi5T/3xl9bTc+kK7mm22g+zfoO\nYfQgyn/dTP7H/0YSRcICfRiYGsakwcEsevgKZkxK5bKUMB6Z2penpl1GenYpr32yG1Fs3fOfQqHC\nP7QnlUV72i3LQLuYbn/6SYxBPkREd/zCfOOuT0WS4KuVe/6Qy2xV1NSTllHI2EEXH5bZkJdGpCGM\nPmGpAEwcGoePl5qVa7PPLDNtK4Yeqah825ZCwekSWfyvXWQfq2TO9MsI0hUDAn5BydzW63pSghN5\nXbkH+8QhnPz8S3567EE2bvses92CXTw3PeLQ3hHMuq0faRlFrFiT2WpdxrA+2K1VWGry23B2zcfj\nprPWO8g+UNzsbM2eRm/QMfnOfuTllLEz7VhHy3E7q3/LQ61SMG5I7AW3K6otZfvJPYyMG9L4d/HS\nqrhjfHd+2XWC9OyG4RVbeQWmzCyChw9rs7a3PtvHjoPFPDltAP27h1JTloWPfxdUGh9UShWPD72f\n7sGJLPXP49uRAShKK1Et+IB/PfVnXt2/jOPVBeeUeXmvCP50dQpfbDjC2m2tM43emIhK40tl8d62\nnmKz8HhHSvrWfFwukb4DPTMbuDUkJIWQ2ieStA25DBgS65G5fB1BndXBD1uPMW5wLL4Xibz5d8ZX\n+OkMjEu88oz3xw2OZUtGIUs+3cubs0dSs2ULglJJ4OBBbdL2/eY8ft55nIdv7sOQnhGILju1FTmE\nxY9u3Ean1jF3xMNYnTaUggKFUyTzi//Q4/PVaJQmXtK9Tpg+hIFRfbkm6b/73TgykcJyC+98sY/w\nIG96JQa3SJsgKDCG9aGyaA9R3a5BofCsLTx6tTkdLrb9mkev/lEY/C/ci9beDB2VSG2NlewDxR0t\nxW2s3ZaPze7i2isunGMkuyyXHSf3ckvPa9GelUVZoRB45Oa+WG1Onl6axsn1mwjo17dVM8RFUUKS\nJD5bn8PSr/Yz5rIuXDWwCwDVpQcRRQcBob3O2U+n0qJWqlFqtfS8bRqJf76XpJxarllbSsLaLNZs\n+ISNuWk4LRagoRfygcm9SI4z8sq/06mqbfm8yKCowTgddZSf3N7ifVuKR023b9cJLGYbl1/pmUQz\nbSE0wkBUTADpW9unHe9prDYn3/yayxV9IwkOaPoGZ3PaWb77P8T4RzE85vy1V4jRmxfuuxyduQpn\n/jEO+8c3u5PB5RL5bH0Of17wMzfN+Y6b5n7PijVZ3DymG49M7dPYlK0o2o2PXww6n4vXSmHjx+J7\nx63E6sNJLoVbf6rB8tzr7Jg+g7JNDWNsKqWCWbf1RxRh0ce7Wtwp4uUbSkBYb4qP/oLoan16/ebg\nMdOJLpG0Dbl07xFGUKjeU4dpE/0vj+HYkXKKC1qW4rwz8u+12Zgsdm65qumZFDannSXbllNQW8z9\nl91xwRne3boE8ECCFVGl5sMjCl77ZPdFL2RRlHjnyww+/iGLpC4BTB7VletHJPDifUO4c0Jyo+FE\n0UltZS7+IanNPj9Nj1R6LVpA/3ffJvS6qykL0lKdEELOq/+gaM2PAAT5e/HXW/tyILeCPa2Is41I\nuAqH3Uzp8S0t3rcleMx0WfuLqaqoY+go9yx35Al69I3EL8CL337uXGtSt5QTJbV8+2suU6/qRmTw\n+ZuBlXXVPPnTy+wuOsAjg2c0GWN5GqfZTPH3a4gYO5oHbxvIhvSTrP6t6bhISZJY/t1Bftqez4NT\n+jDr9v7cNq47d4xPpk+3M/PTWGqOI4mOVuW3VGq1dJv+J4z33cGKfi4MY0aQ9+575C79J6LDQb+k\nEKJD9fy0veUtGJ1PCEERl1F8bAOuNqTuuBgee2LM2HWCyC7+7TqFp6UolQqGje7K959nUFpkIiTc\n0NGSWsXHP2QR6O/FDSOavog/2vs5FnsdC8fOIdrv4qF4BV99g+RwEn3TFBKMARw5Uc373x7g0PEq\neiUGY3e40GmURIXosVgd/LQ9n637i5gxKZVxgy9s6NrKXBQqXZtylIxNHMHq7J/57TI/bki8l6PL\nllO5YxdxM6YzdlAMH353kKpaKwF6XYvKDY8fTUXhLkqPbyH8d5087sQjNZ213kFuThkpfTomzrIl\n9LksGr8AL75ZtfeSDIbOLzKxdX8RU8ckoTnPEtEAB0tz2Hoindt6Xd8sw9mrayj8bg3hV09AY2yY\na3b3dT2YNjGF/CITb362l2XfHuCNT/fyxJu/8fyybew/Us6T0wZww5UXr71qK3PR+8e1acKqTqXl\nuuSxbDq+HYb1o89ri/GJj+PQ/71Gf6kIQRD4eUfLMzlrvAIIihpIybFNuBz1F9+hFXikpsvKKEJ0\nSaT0CvdE8W5FqVJw0/QBfPR2Gl+t3MPN0wdcUnGZn67PIchPx6gmFtG0uxx8kL6KBGNMYyqFi3Hy\n8y8bEhzd+N/EtyqlgsmjujJ5VFfqbU60aiVWu5Piijp8vNQE+ematQKu6HJgqcknMnF8807wAlyV\nMJxvs9fx2YHveGTIDJLnPkn2wv/j5FtvMW7E7axJO8aNVya2eGXesLhRlBfspCT/NyISx7ZZ59m4\nvaaTJImdW44S3y0Yv4BzA2c7IxHR/txwW18OHSzm60/24LpEIlUKysxs3lvAjSO7om5irPGHnA0U\n1Jbw5wG3N6ZHuBCmrGyKf/iRiOsmoTacv7ntpVWhUAh469TER/oRavRu9oVtrs5HEp1uWa9Aq9Iw\nJXUim4/vJLssF0GhoNtjM/GKjqLPnu+oraghLaOoxeVqdP4ERw+h+NhGrHXlbdZ5Nm433YmjlRQX\nmBh4hWcWIfQU3XuGc8OtfTm4t5CvV7ZfHF5b+HhNFn6+WsY28Qxltlv4Onsto+KHEhdw8eAEp9nC\nocWv4puYQNSUG90tF2hYZVWp8sJL755W0Jj4K4jzj2b57v/gEl0otVqSZj8K9RamVf3KdxtatzRW\nRMJY1Bpfjh1YhSS63KL1NG41nSRKrPsui6BQX7q6admr9qRn/yiuP2W8A3vODTnqTGQerWBLRiHT\nr05Be55nub1FmTy/4R+IosiUlInNKjNv2Qe46uro9thfz7uqqjswlR9Cb0xwWxZnhULBXf2mcqzm\nJAt+fYvcyny8IiJIeeZpAuvKGbjtPxxI29ficpUqHXE9b8VSc4KCw2vcorVRszsL2739OAX5VUy8\nsecl9Vz0e3r0iySpRxjrVmdis3bOjhVJklixJovYcAMj+0ef89mnB1bz8q9L0CjVPD3iYYzeF+9B\nth/IpGzDRuLuuQtdqGdumHWmk9SZTmIM7+fWcrsHJzDr8j9TVlfBs+v/j2NVJzAkd6fngpfwlezU\nLHyBzJdeboxgaS6+AXFEdp1ASf6vVJXsd5tet5murLiWdasP0mtAFLEdNFHVXYy9NhVrnYONa9t3\nGn9zSc+h61jdAAAQEklEQVQu5WBeBXdOTEbxu5ubKIos3/0pnx9cw9Qek3hp9ON0C4q/QEkN2Ktr\nsHz1DQGX9SdktOcyLBcf3YBaa8A/OMXtZQ+M6sPisU8Trg9l4eZ3OFFTiCExAdtf5vJd2DBqDmSR\n8eRcbGUte0YLjRmBf0gq+Qc/xWlvmWmbos2mE0WJzH2FfPL+Dvz8vZhwQ/vnEXQ3AYHejBiXxPZf\n89j0U06ner6rrbPzzhf7SIkzclly6BmfLd35L9bmbuLufrcwOXVis2Z1SJJE3tJ3QZJIfPABj80E\nMVXkUFWSQWTiBI/lttSoNDw1/C94q3Q8/fMi9hVnMmZIPPkh3dk38g5c9VYyX5yPs67u4oWdQhAE\nuiRPQZIkivJ+dovONpnOWu9g+Ztb+HxFOj6+GqbOGNih68C5k8tHJjBibDc2rT3ELz9kX3yHdkCS\nJN74zx4sVieP3trvDIPsLNjHxmNbuW/A7YzrOqLZZZZt+o2Krdvxvm4SmgDP5H90Oa3kH/wM34B4\njBHubVqeTZC3kZfGPEG3wHheTXuPetHMrWOTWHPYiuLO+7CVlpH98kJER/PjK9VaX8LirqT0RBpW\nS9vTOLbadHUWO/96dxsVpWam/+Vy7p55BcYg968T3VEIgsCIcUmMvjqZLeuPsHn9YaRWzk52F7sP\nlbLtQDEPTu5NWOB/v2uzzcJ7u1bSNzyVkXHNX6TRkn+c3LeXEnTFUDS9enhCMk67hSN7luO0W4hN\nvald1oTzUuv46+V3o1GoWZb+CRMuj6VnQhAvrTmJ6Ya7MGVlc2z5Ry0qMzRmOGqtgUM73mpzULTy\nueeee64lOxQVFVFVIvHlv3ZjNlm5477BRMd6duWd5lJeXk5wcMvmUl2M6NgAXE6RX3/KobLcQmJy\nSIsHW92hTxQlFq7YSXigD3dfm/q74GGR17d9QFFtKXNGPHROKvSmcJhMHPzb86gNBpKfforKmhq3\nf3em8kMc3r0Mh81EYt+78Da0Pl1HS787jVJDiE8gX2b+QLCPkemjBlNVa2fl9lJiYkNQbPoRTWAg\nvgkXf+YFEBRKAkJ74XSYGzpWivZgrSvFXJOP2eZDRETzo69a1Rb84l+7ie8WzIQbexDYRIDtHwVB\nEBh9dTKh4Qa++c9eSotMTLixJzEJge2qY0P6CY4Wmlj40LBGwzldTt7Z+THphft5fNj9BHlf/OYn\nOp1UbtvOiU8/x2mpo9eil1F6uX+uo6X6OEf2foivfywxqTeh9Wr/G/Pg6H6MjLuc99I/IcjbyINT\nepMY5cfSLzK4MTwV3nkXTYA/xsual8tTo/MnJmUKIdHDOHn4e2or83A6LKC/pkW6WmW6GY8MIyqm\nffO/dzQ9+kUSHKbnm1V7+OjtNJJSQ5lwY892mZybebSCt7/IYFjvCFLiGsx+tOoEH+/9gqzyIzw8\n+E8MiDx3MuhpXPX1VGzfifnIESrStmKvqESf3J3UmQ/jFe7+UL3aylxy936Etz6CxL4zUCg7bv2I\ne/vfSmV9FQs3v8NDg6YzbnB/4iP9WPCBGkWdGRa+Qs+XnsPQvfmLy3jpw+ja7+7G39PT01ukqVWm\n+18z3GlCIwzc+9fhHNhbwE/fZvL2oo2k9A4nKERPt5QQj8wbzDlexQvLttE12p+Hp/Zmx8m9rMn5\nhcyywwR6BTDnigfpFXZuAljJ5cJaUkL5b1soXP09ztpatCEh+PftS/jE8c1uVrUESRKpLNpNfuYX\n+PrFEN9nWocaDmjIvTLsAd7Y9gGvpS0jLSqdm1KvZtHMEbz8TxXq9C8Q5z5L8Mgr8e+WgD4pCZ9Y\nz64o1SGLQnqK9lx0sb7Ozq/rDnPsSDkVpWacTpHAYB9iE4OIjg0gNjHonFrwtD5JkkBqWJ5ZdDgQ\nFAoE5Znd6JIk8f22wyzfsAm/MDO9kvzJLc3BVl5Oom8kQyP60FUTimiqxVZaRt3xE9irq3HU1CBa\nbTjNZiSXC4VGQ/DIEURNvgFd6JlDDOfT1lrM1fmcyP6aOtNJAkJ7EdvjFrcarq36JEnit/wd/Gvf\nl1RbTfhp9fQI7Y6iLATxq+0kWk7iZ6tFQMK/T2/CJ12Nf6+eKDSai5bdUk+4xXSSKCK5XCAICEpl\n4zOHJEkgik28Sg0LN0hSQ6+gJJ56Pf2zCAgggEMSMVlNVFuqqbXWYrGZKa2toMpSic1hxeaw43Ta\nsVttSArQarSEGUKJ8Asn0hCGES8ULheC04VgdyI6nEh2O6LdhsvuQHI4kZyOBgMIApJSiaPWhLOi\nCmetGcUpUwiCAEoFKBSgVCAJCiRBQPAPoNjpS4lZS1m9Fouj4WLzcZkwSGZcSi0upRpJkJBQYXFp\nsKNGkERUoh1fWyVhdUcJthcgqpTUKyScogOF5EQhSaglBSoXaJuIkFEH+OMTE4PGGIDKYEDp5YVa\n74suPBzfromo9RevgZu6qCVJxFx1DNFlaxgcFgRMFYex1ZUjSa5TY24ClupjePmGE939OvRG96fn\ncNcN1e60k1GSTXb5EXYW7KOothSNUoPeFUnNUV+6HLMwqi4Pn5pSFDodQUOHEDZ+XGPLQLTbER1O\noME2Sp2OPfv3e950ppf/jkKSUIgSggRnD6dKnPteZ8ShBKdSwKUUcCnApWxQrXRJWLUKzF4K6nUN\nPZUKERRSw/kqRBpeJQmFCAaLC4UIktDwz67QctIYRZ0qHEnUo5RsqFxOFKICBAeCog5RacWuVmBV\naREcwQjOEMCGSqhCJZpRCwJ6nZq4IB8i/DR4e6nQhQShDQ5GqdMhqFRo/P1RGfRtSnUuSSKS6CQ7\n6yCJiXGILjtORx1Ou4V6czEVhbuwnRVpr/UOxscvGkGhQhKdSKIT34B4gqMGe2zg2xOtGEmSOFFT\nyM6CfewqyCC3Kh8BAaHegP9JNQnFZnqWlKC32xARUHCuVQSVCu3cJzxvumO7VzYc/tRFJtHw2lAg\nCJKE1FBJIZ2qraTTNjxdC54+8UZ3Cg3/hIZXAUXDqyQ0fBEoEVCBpEQQFA3jPYKASqlEqWp4tdnt\naLRaBKmhxhQQcQkOnEoRSRBwChIi4GrUdfrgDc29ht0kJEApKBCEht9FRCSkBhUKAaVCiVIhoFYp\nUSlO/SkkAVGUcLpcCChQCAokRASh4RYkCGCz2dDpNKfKVqAUFKiUCkRJQuEScNWJiA4Ju82JJEo4\nHK7GzMVKpYBSqUCpVKBQKlAqBRQKAaVKgUIh/DfWVTr9n4SAE0FwAq5TP5/5iuBCITQdQS9JKpxE\n4yAZEQMSulPfnrYll4xbKC8rI8jNQxpnU++op6i2jFJLBeWWSmwuG4gieouIziYiCSAqhMabK4DW\nLjFhyhjPm05GRuZMPGo6GRmZtvHHSG0sI3MJIZtORqadkU0nI9POyKaTkWlnLjjAs2fPHlatWoWv\nry9GoxGdTkdBQQFms5k5c+YQ4KH5V80lPz+f119/HaPRSGpqKpWVlZ1KH8Ds2bMZNWoURUVFnUpb\nYWEhf/nLX0hOTiYkJASDwdCp9BUUFPD222/j6+uLwWDoVNfev//9bw4cOIDD4WD37t3cfvvtLdJ2\nwak9OTk53HbbbYwZM4Z3332XkpISFi9ejMFgYP369QwY0LzobE9RUFDAqFGjGD9+PK+++mqn0/fh\nhx9isViIiopi48aNnUrbL7/8QkFBAd7e3gwdOpR169Z1Kn1vvPEGERERlJSUMHToUH766adOo69X\nr16MGTOG/fv388ADD7B69eoWabtgTTdiRMMM5KVLlzJp0iR27doFQFhYGGVlbZ9B21ZSU1MpKSnh\nvvvuY9CgQRw/3pDRtzPo++WXX9Dr9fTp0wdRFAkMDOw02qDhwhk6dChBQUFMnz6d6OiGBEedRV9+\nfj5Tpkzhnnvu4a677qJLl4YltjqLvry8PJxOJxERES3+217wmc5isfD000/Tp08fJk2aRHV1NQDF\nxcWEhHR8ir3MzEy0Wi3vv/8+Bw4coKqqCugc+lavXs3+/fv56quv+PTTT6moqOg02qAhrMputwOg\n0+koLS0FOo++4OBgfHx8UCqV6HS6TnftrVy5kjvvvJPAwMAWa7vg4PjcuXM5fvw4ERERKJVKUlJS\nyMvLw2Qy8fzzz+PbioUC3UlGRgbLli0jLCwMtVpNWFhYp9IH8PXXX6PRaKioqOhU2g4ePMh7771H\nUFAQSUlJWK3WTqUvNzeXJUuWoNfrGTRoEFVVVZ1K31133cXy5csB+Pjjj1ukTY5IkZFpZ+QhAxmZ\ndkY2nYxMOyObTkamnZFNJyPTzsimk5FpZ2TTuYmCggIeeeQRAMaOHcu0adOYNm0at9xyC+vXrz/j\n84tx/fXXX3yj3/H1118D8Oabb7Jp06Zm7/fUU09x0003UVFRwX333cf06dMxm80APPPMM43jePPm\nzWsMlJBpO3+MhQc6CacTMhkMBlasWAE0BBhMmzaNJUuWNHtxjpYu4vHxxx+32Kinj7NgwQKKiooY\nNmwYkZGRpKWlERkZSdeuXdGcyoT10ksvMWXKlBaXL3N+ZNN5gN8PfdbU1ODt7X3RlX8+++wzVq1a\nRVxcHE5nQ9av9evX895776FQKJg5cyaDBg1i0qRJREdHU1xczGOPPYbFYuHo0aO89tprqNVqPvnk\nE95//330ej1vvPEGa9euZcWKFYiiyEMPPcTw4cPPObaXlxc2m436+noMBgOrVq3imWeece+XItOI\n3Lz0AGazubF5+cILLzB37tyL7rNy5Uo+++wzHnvsMex2O5Ik8c477/Dxxx/zwQcfsGTJEgBKS0tZ\ntGgRy5cvZ+nSpYwbN474+HgeffRRAPr168eKFSswGAwcPHiQH3/8kXnz5rFs2bJGM59NQkIC1dXV\n7NmzBz8/PxISEvjkk0+YP39+Y3NTxn3INZ0H0Ov1jc3L0xQUNL2ccnl5OWFhYSgUCiIiIggICKCy\nspKTJ09y9913I0kSVVVV2O124uLiGsOMbDYbcGbNmpTUkB7c398fm83G7NmzeeuttygrK+Pmm29u\nUsPs2bOBhme5mTNn8t5773HDDTewevVqbr311k61Rt+ljmw6D9DUBdrU+/7+/hQXF+NwOKisrKSq\nqgqj0Uh8fDwfffQRTqeT5cuXo9FoyM/Px2q1UldXh/JUVuimngElSeLbb7/lmWeeQaPRMGXKFMaP\nH9+k7szMTOLj41GdyqMpCEKjsWXch2w6D9CUCbZv397YIZGcnMyLL74IgEqlYsaMGUydOpXY2Fj8\n/PwQBIE777yT22+/HavVyi233NK47cyZM6msrOTJJ58EICYmhvnz5+Pn53eGBkEQSEpK4o477sDX\n15c777zzgrpXrVrF3LlzG6P6X3rpJRYvXnzBc5JpOXLA8yXGpEmTWL16tVvKmjNnDjNmzKBr164X\n3Xby5Ml88cUXbjnu/zpyR8r/OHPnzm2ch9gU8+bNa5wPKNN25JpORqadkWs6GZl2RjadjEw7I5tO\nRqadkU0nI9POyKaTkWlnZNPJyLQz/w/QoRAVwlSiMAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbc8bd9f2b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(9/2.54,4/2.54))\n",
    "sns.set(font='Arial',style='whitegrid',context='paper')\n",
    "plt.plot(np.linspace(-17.423,119.651,3576),np.sum(cbc5_volume_density[:,cbc5_prediction==1],axis=1)/np.sum(cbc5_volume_density[:,cbc5_prediction==1]),label='CBC 5T')\n",
    "plt.plot(np.linspace(-17.423,119.651,3576),np.sum(cbc5_volume_density[:,cbc5_prediction==2],axis=1)/np.sum(cbc5_volume_density[:,cbc5_prediction==2]),label='CBC 5O')\n",
    "plt.plot(np.linspace(-17.423,119.651,3576),np.sum(cbc5_volume_density[:,cbc5_prediction==0],axis=1)/np.sum(cbc5_volume_density[:,cbc5_prediction==0]),label='CBC 5I')\n",
    "plt.plot(np.linspace(-17.423,119.651,3576),np.sum(off_sac_volume_density,axis=1)/np.sum(off_sac_volume_density[700:,:]),label='OFF SAC')\n",
    "plt.plot(np.linspace(-17.423,119.651,3576),np.sum(on_sac_volume_density,axis=1)/np.sum(on_sac_volume_density[:2300,:]),label='ON SAC')\n",
    "plt.xlabel('IPL depths [%]')\n",
    "plt.yticks([])\n",
    "plt.xlim(20,70)\n",
    "plt.legend(bbox_to_anchor=(0.6, 1.05))\n",
    "# plt.savefig('figures/type5_profiles.svg',bbox_inches='tight',dpi=300,transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot dendritic and axonal mosaics of the new CBC5 subtypes (Supp. figure 7C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#5T\n",
    "plot_mosaic(CBC5T-390,'IPL')\n",
    "plot_mosaic(CBC5T-390,'OPL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#5O\n",
    "plot_mosaic(np.where(cbc5_prediction==2)[0]+144,'IPL')\n",
    "plot_mosaic(CBC5O-390,'OPL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#5I\n",
    "plot_mosaic(CBC5T-390,'IPL')\n",
    "plot_mosaic(CBC5I-390,'OPL')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def calculate_overlap(selection,layer):\n",
    "    total_nodes=np.empty((0,2))\n",
    "    polygons=[]\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",
    "            if layer=='OPL':\n",
    "                nodes=(nodes*[16.5,16.5,25]+[0,3250,15000]).astype(int)\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,19000]).astype(int)\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([],[])\n",
    "        if nodes_cell.shape[0]>2:\n",
    "            nodes_cell=nodes_cell[:,1:]/[25,25]\n",
    "            total_nodes=np.concatenate((total_nodes,nodes_cell),axis=0)\n",
    "            hull=ConvexHull(nodes_cell)\n",
    "            polygons.append(Polygon(nodes_cell[hull.vertices,:]))\n",
    "    hull=ConvexHull(total_nodes)\n",
    "    overlaps=[]\n",
    "    area=[]\n",
    "    for i in range(len(polygons)):\n",
    "        area.append(polygons[i].area)\n",
    "        for j in range(i+1,len(polygons)):\n",
    "            overlaps.append([BC_IDs[selection,0][i],BC_IDs[selection,0][j],polygons[i].intersection(polygons[j]).area])\n",
    "    overlaps=np.array(overlaps)\n",
    "    area=np.array(area)\n",
    "    return(area,overlaps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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,19000]).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",
    "    plt.figure(figsize=(15,15))\n",
    "    plt.imshow(im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}