docs: def plot_population(ax,
                      params,
                      aspect='tight',
                      isometricangle=0,
                      plot_somas = True, plot_morphos = False,
                      num_unitsE = 1, num_unitsI=1,
                      clip_dendrites=False,
                      main_pops=True, 
                      Y = None,
                      big=True,
                      title='cell positions',
                      rasterized=True):
      '''
      Plot the geometry of the column model, optionally with somatic locations
      and optionally with reconstructed neurons
      
      kwargs:
      ::
          ax : matplotlib.axes.AxesSubplot
          aspect : str
              matplotlib.axis argument
          isometricangle : float
              pseudo-3d view angle
          plot_somas : bool
              plot soma locations
          plot_morphos : bool
              plot full morphologies
          num_unitsE : int
              number of excitatory morphos plotted per population
          num_unitsI : int
              number of inhibitory morphos plotted per population
          clip_dendrites : bool
              draw dendrites outside of axis
          mainpops : bool
              if True, plot only main pops, e.g. b23 and nb23 as L23I
          Y : None, or string
              if not None, plot only soma locations of Y
          big : bool
              if False: leave out labels and reduce marker size
      
      return:
      ::
          axis : list
              the plt.axis() corresponding to input aspect
      '''
  
      name_to_id_mapping={'L6E':(3,0),
                          'L6I':(3,1),
                          'L5E':(2,0),
                          'L5I':(2,1),
                          'L4E':(1,0),
                          'L4I':(1,1),
                          'L23E':(0,0),
                          'L23I':(0,1)
      } 
      
      
      # DRAW OUTLINE OF POPULATIONS 
      
      ax.xaxis.set_ticks([])
      ax.yaxis.set_ticks([])
  
      #contact points
      if big:
          ax.plot(params.electrodeParams['x'],
              params.electrodeParams['z'],
              '.', marker='o', markersize=2, color='k', zorder=0)
      else: 
          ax.plot(params.electrodeParams['x'],
              params.electrodeParams['z'],
              '.', marker='o', markersize=0.5, color='k', zorder=0)
  
      #outline of electrode       
      x_0 = params.electrodeParams['r_z'][1, 1:-1]
      z_0 = params.electrodeParams['r_z'][0, 1:-1]
      x = np.r_[x_0[-1], x_0[::-1], -x_0[1:], -x_0[-1]]
      z = np.r_[100, z_0[::-1], z_0[1:], 100]
      ax.fill(x, z, fc='w', lw=0.1, ec='k', zorder=-0.1, clip_on=False)
  
      #outline of populations:
      #fetch the population radius from some population
      r = params.populationParams['p23']['radius']
  
      theta0 = np.linspace(0, np.pi, 20)
      theta1 = np.linspace(np.pi, 2*np.pi, 20)
      
      zpos = np.r_[params.layerBoundaries[:, 0],
                   params.layerBoundaries[-1, 1]]
      
      layers = ['L1', 'L2/3', 'L4', 'L5', 'L6']
      for i, z in enumerate(params.layerBoundaries.mean(axis=1)):
          if big:
              ax.text(r, z, ' %s' % layers[i], va='center', ha='left')
  
      for i, zval in enumerate(zpos):
          if i == 0:
              ax.plot(r*np.cos(theta0),
                      r*np.sin(theta0)*np.sin(isometricangle)+zval,
                      color='k', zorder=-r, clip_on=False)
              ax.plot(r*np.cos(theta1),
                      r*np.sin(theta1)*np.sin(isometricangle)+zval,
                      color='k', zorder=r, clip_on=False)
          else:
              ax.plot(r*np.cos(theta0),
                      r*np.sin(theta0)*np.sin(isometricangle)+zval,
                      color='gray', zorder=-r, clip_on=False)
              ax.plot(r*np.cos(theta1),
                      r*np.sin(theta1)*np.sin(isometricangle)+zval,
                      color='k', zorder=r, clip_on=False)
      
      ax.plot([-r, -r], [zpos[0], zpos[-1]], 'k', zorder=0, clip_on=False)
      ax.plot([r, r], [zpos[0], zpos[-1]], 'k', zorder=0, clip_on=False)
      
      if big:
          #plot a horizontal radius scalebar
          ax.plot([0, r], [z_0.min()]*2, 'k', lw=1, zorder=0, clip_on=False)
          ax.text(r / 2., z_0.min()-100, '$r$ = %i $\mu$m' % int(r), ha='center')
      
          #plot a vertical depth scalebar
          ax.plot([-r]*2, [z_0.min()+50, z_0.min()-50],
                  'k', lw=1, zorder=0, clip_on=False)
          ax.text(-r, z_0.min(), r'100 $\mu$m', va='center', ha='right')
      
      ax.set_yticks([])
      ax.set_yticklabels([])
      
      #fake ticks:
      if big:
          for pos in zpos:
              ax.text(-r, pos, '$z$=%i-' % int(pos), ha='right', va='center')
          ax.set_title(title, va='bottom')
     
      axis = ax.axis(ax.axis(aspect))
  
  
      def plot_pop_scatter(somapos, marker, colors, i):
          #scatter plot setting appropriate zorder for each datapoint by binning
          pitch = 100
          for lower in np.arange(-600, 601, pitch):
              upper = lower + pitch
              inds = (somapos[:, 1] >= lower) & (somapos[:, 1] < upper)
              if np.any(inds):
                  if big:
                      ax.scatter(somapos[inds, 0],
                                 somapos[inds, 2] - somapos[inds, 1] * np.sin(isometricangle),
                                 s=10, facecolors=colors[i], edgecolors='gray', linewidths=0.1, zorder=lower,
                                 marker = marker, clip_on=False, rasterized=rasterized)
                  else:
                      ax.scatter(somapos[inds, 0],
                                 somapos[inds, 2] - somapos[inds, 1] * np.sin(isometricangle),
                                 s=3, facecolors=colors[i], edgecolors='gray', linewidths=0.1, zorder=lower,
                                 marker = marker, clip_on=False, rasterized=rasterized)
                      
  
  
  
      # DRAW UNITS
      pop = zip(*params.mapping_Yy)[0]
      
      #plot a symbol in each location with a unit 
      if plot_somas:
          if main_pops:
              colors = phlp.get_colors(np.unique(pop).size)
  
              #restructure 
              E, I = zip(*params.y_in_Y)
              pops_ = []
  
              if Y is None:
                  for i in xrange(len(E)):
                      pops_.append(E[i])
                      pops_.append(I[i])
              else:
                  ids = name_to_id_mapping[Y]
                  if ids[1] == 0:
                      pops_.append(E[ids[0]])
                  if ids[1] == 1:
                      pops_.append(I[ids[0]])
  
              for i, pops in enumerate(pops_):
                  layer = np.unique(pop)[i]
  
                  if layer.rfind('p') >= 0 or layer.rfind('E') >= 0:
                      marker = '^'
                  elif layer.rfind('b') >= 0 or layer.rfind('I') >= 0:
                      marker = '*'
                  elif layer.rfind('ss') >= 0:
                      marker = 'o'
                  else:
                      raise Exception
                  
                  #get the somapos
                  somapos = []
                  for j, lname in enumerate(pops):
                      fname = glob.glob(os.path.join(params.populations_path, '%s*somapos.gdf' % lname))[0]
                      if j == 0:
                          somapos = np.loadtxt(fname).reshape((-1, 3))
                      else:
                          somapos = np.r_['0, 2', somapos, np.loadtxt(fname).reshape((-1, 3))]
                  somapos = somapos[::5, :]
                  if Y is None:
                      plot_pop_scatter(somapos, marker, colors, i)
                  else:
                      plot_pop_scatter(somapos, marker, colors, ids[0]*2+ids[1])
              
          else:
              colors = phlp.get_colors(len(pop))
              i = 0
              for layer, _, _, _ in params.y_zip_list:
                  
                  #assign symbol
                  if layer.rfind('p') >= 0 or layer.rfind('E') >= 0:
                      marker = '^'
                  elif layer.rfind('b') >= 0 or layer.rfind('I') >= 0:
                      marker = '*'
                  elif layer.rfind('ss') >= 0:
                      marker = 'x'
                  else:
                      raise Exception
                      
                  #get the somapos
                  fname = glob.glob(os.path.join(params.populations_path, '%s*somapos.gdf' % layer))[0]
                  somapos = np.loadtxt(fname).reshape((-1, 3))
              
                  plot_pop_scatter(somapos, marker, colors, i)
                  
                  i += 1
      
  
      #plot morphologies in their appropriate locations
      if plot_morphos:
          if main_pops:
              colors = phlp.get_colors(np.unique(pop).size)
  
              #restructure 
              E, I = zip(*params.y_in_Y)
              pops_ = []
              for i in xrange(len(E)):
                  pops_.append(E[i])
                  pops_.append(I[i])
  
              for i, pops in enumerate(pops_):
                  layer = np.unique(pop)[i]
                  
                  #get the somapos and morphos
                  somapos = []
                  for j, lname in enumerate(pops):
                      fname = glob.glob(os.path.join(params.populations_path, '%s*somapos.gdf' % lname))[0]
                      if j == 0:
                          somapos = np.loadtxt(fname).reshape((-1, 3))
                      else:
                          somapos = np.r_['0, 2', somapos, np.loadtxt(fname).reshape((-1, 3))]
  
  
                  #add num_units morphologies per population with a random z-rotation
                  if layer.rfind('p') >= 0 or layer.rfind('ss') >= 0 or layer.rfind('E') >= 0:
                      num_units = num_unitsE
                  else:
                      num_units = num_unitsI            
                  
                  if num_units > somapos.shape[0]:
                      n = somapos.shape[0]
                  else:
                      n = num_units
  
                  #find some morphos for this population:
                  morphos = []
                  for fname in params.m_y:
                      if fname.rfind(layer) >= 0:
                          morphos.append(fname)
                                  
                  #plot some units
                  for j in range(n):
                      cell = LFPy.Cell(morphology=os.path.join(params.PATH_m_y,
                                              np.random.permutation(morphos)[0]),
                                       nsegs_method='lambda_f',
                                       lambda_f=10,
                                       extracellular=False
                                      )
                      cell.set_pos(somapos[j, 0], somapos[j, 1], somapos[j, 2])
                      cell.set_rotation(z=np.random.rand()*np.pi*2)
      
                      #set up a polycollection
                      zips = []
                      for x, z in cell.get_idx_polygons():
                          zips.append(zip(x, z-somapos[j, 1] * np.sin(isometricangle)))
                      
                      polycol = PolyCollection(zips,
                                               edgecolors=colors[i],
                                               facecolors=colors[i],
                                               linewidths=(0.5),
                                               zorder=somapos[j, 1],
                                               clip_on=clip_dendrites,
                                               rasterized=rasterized)
      
                      ax.add_collection(polycol)
                  
                  i += 1
  
  
          else:
              colors = phlp.get_colors(len(pop))
              i = 0
              for layer, morpho, depth, size in params.y_zip_list:
      
      
                  #get the somapos
                  fname = glob.glob(os.path.join(params.populations_path, '%s*somapos.gdf' % layer))[0]
                  somapos = np.loadtxt(fname).reshape((-1, 3))
      
                  
                  #add num_units morphologies per population with a random z-rotation
                  if layer.rfind('p') >= 0 or layer.rfind('ss') >= 0 or layer.rfind('E') >= 0:
                      num_units = num_unitsE
                  else:
                      num_units = num_unitsI            
                  
                  if num_units > somapos.shape[0]:
                      n = somapos.shape[0]
                  else:
                      n = num_units
                  
                  #plot some units
                  for j in range(n):
                      cell = LFPy.Cell(morphology=os.path.join(params.PATH_m_y, morpho),
                                       nsegs_method='lambda_f',
                                       lambda_f=10,
                                       extracellular=False
                                      )
                      cell.set_pos(somapos[j, 0], somapos[j, 1], somapos[j, 2])
                      cell.set_rotation(z=np.random.rand()*np.pi*2)
      
                      #set up a polycollection
                      zips = []
                      for x, z in cell.get_idx_polygons():
                          zips.append(zip(x, z-somapos[j, 1] * np.sin(isometricangle)))
                      
                      polycol = PolyCollection(zips,
                                               edgecolors=colors[i],
                                               facecolors=colors[i],
                                               linewidths=(0.5),
                                               zorder=somapos[j, 1],
                                               clip_on=clip_dendrites,
                                               rasterized=rasterized)
      
                      ax.add_collection(polycol)
                  
                  i += 1
      
      return axis