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
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