def stellar_massftn():
    gadget2msun=10.e10
    boxsize = 47.0
    max_mag=-16.
    min_mag = -23.
    nbins=14
    hubble_h = 0.7
    #subfind_folder = "/mnt/lustre/scratch/cs390/nIFTy/62.5_dmSF/outputs/"
    #ahf_folder = "/mnt/lustre/scratch/cs390/nIFTy/62.5_dm/outputs/"
    
    firstfile = first1
    lastfile = last1

    filter = LGalaxyStruct.properties_used
    filter['DiskMass'] = True
    filter['BulgeMass'] = True

    #file_prefix = "SA_z0.00"    
    (nTrees,nGals,nTreeGals,gal) = read_lgal.readsnap_lgal(folder1,file_prefix,first1,last1,filter)
    massf = gadget2msun*gal['DiskMass']+gadget2msun*gal['BulgeMass']
    mass = numpy.log10(massf)
    stellarmass = pylab.histogram(mass,bins=20,range=(9.0,14.0))
    print stellarmass
    massftn_y = stellarmass[0]
    massftn_x = []
    for i in range(len(stellarmass[0])):
        massftn_x.append((stellarmass[1][i]+stellarmass[1][i+1])/2.)

    delta_logM = massftn_x[1]-massftn_x[0]
    pylab.rc('text', usetex=True)
    fig = pylab.figure()
    ax = fig.add_subplot(111)
    ax.plot(massftn_x,massftn_y/boxsize1**3./delta_logM,'r-',label=label1)

    firstfile = first2
    lastfile = last2

     
    (nTrees,nGals,nTreeGals,gal) = read_lgal.readsnap_lgal(folder2,file_prefix,first2,last2,filter)
    massf = gadget2msun*gal['DiskMass']+gadget2msun*gal['BulgeMass']
    mass = numpy.log10(massf)
    stellarmass = pylab.histogram(mass,bins=20,range=(9.0,14.0))
    print stellarmass
    massftn_y = stellarmass[0]
    massftn_x = []
    for i in range(len(stellarmass[0])):
        massftn_x.append((stellarmass[1][i]+stellarmass[1][i+1])/2.)
    ax.set_xlabel(r"$\log(M_\star/M_\odot$ $h)$")
    ax.set_ylabel(r"galaxies$/(Mpc^3 h^{-3})/\Delta \log(M_\star/M_\odot$ $h)$")
    ax.plot(massftn_x,massftn_y/boxsize2**3./delta_logM,'b-',label=label2)

    print "Stellar mass"
    for i in range(len(massftn_x)):
        print massftn_x[i],"\t",massftn_y[i]/boxsize**3./delta_logM
    ax.set_yscale("log")
    ax.legend(loc='upper right',ncol=1, fancybox=True)


    #pylab.show()
    pylab.savefig('stellar_mass.pdf',bbox_inches='tight')

    def get_bcc_pz(self,filename_lenscat):

        if self.prob_z == None:

            # filename_lenscat = os.environ['HOME'] + '/data/BCC/bcc_a1.0b/aardvark_v1.0/lenscats/s2n10cats/aardvarkv1.0_des_lenscat_s2n10.351.fit'
            # filename_lenscat = os.environ['HOME'] + '/data/BCC/bcc_a1.0b/aardvark_v1.0/lenscats/s2n10cats/aardvarkv1.0_des_lenscat_s2n10.351.fit'

            if 'fits' in filename_lenscat:
                lenscat = tabletools.loadTable(filename_lenscat)
                if 'z' in lenscat.dtype.names:
                    self.prob_z , _  = pl.histogram(lenscat['z'],bins=self.grid_z_edges,normed=True)
                elif 'z-phot' in lenscat.dtype.names:
                    self.prob_z , _  = pl.histogram(lenscat['z-phot'],bins=self.grid_z_edges,normed=True)

                if 'e1' in lenscat.dtype.names:

                    select = lenscat['star_flag'] == 0
                    lenscat = lenscat[select]
                    select = lenscat['fitclass'] == 0
                    lenscat = lenscat[select]
                    select = (lenscat['e1'] != 0.0) * (lenscat['e2'] != 0.0)
                    lenscat = lenscat[select]
                    self.sigma_ell = np.std(lenscat['e1']*lenscat['weight'],ddof=1)

            elif 'pp2' in filename_lenscat:

                pickle = tabletools.loadPickle(filename_lenscat,log=0)
                self.prob_z =  pickle['prob_z']
                self.grid_z_centers = pickle['bins_z']
                self.grid_z_edges = plotstools.get_bins_edges(self.grid_z_centers)

def spike_psth(spike_time_ms, t1_ms = -50., t2_ms = 250., bin_ms = 1):
  """."""
  N_trials = len(spike_time_ms)
  t2_ms = pylab.ceil((t2_ms - t1_ms) / bin_ms)*bin_ms + t1_ms
  N_bins = (t2_ms - t1_ms) / bin_ms
  
  spike_count_by_trial = pylab.zeros((N_trials,N_bins),dtype=float)
  if N_trials > 0:
    all_spikes_ms = pylab.array([],dtype=float)
    for trial in range(len(spike_time_ms)):
      if spike_time_ms[trial] is None:
        continue
      idx = pylab.find((spike_time_ms[trial] >= t1_ms) & 
                       (spike_time_ms[trial] <= t2_ms))
      spike_count_by_trial[trial,:], bin_edges = \
        pylab.histogram(spike_time_ms[trial][idx], bins = N_bins, 
                        range = (t1_ms, t2_ms))
      
    spike_rate = 1000*spike_count_by_trial.mean(axis=0)/bin_ms
  else:
    spike_rate = pylab.nan

  dummy, bin_edges = \
    pylab.histogram(None, bins = N_bins, range = (t1_ms, t2_ms))
  bin_center_ms = (bin_edges[1:] + bin_edges[:-1])/2.0

  return spike_rate, spike_count_by_trial, bin_center_ms

def plot_hist(X,Y,title,name):
    # get list of tracks and list of labels
    xs = X.values
    ys = Y.values
    ys = pl.reshape(ys,[ys.shape[0],])
    
    pl.figure(figsize=(15, 6), dpi=100)
    for i in range(many_features):
        if (i==2):
            counts0, bins0 = pl.histogram(xs[ys==0,i],100,range=(0.,0.08))
            counts1, bins1 = pl.histogram(xs[ys==1,i],100,range=(0.,0.08))
        elif (i==5):
            counts0, bins0 = pl.histogram(xs[ys==0,i],100,range=(1,5))
            counts1, bins1 = pl.histogram(xs[ys==1,i],100,range=(1,5))
        elif (i==6):
            counts0, bins0 = pl.histogram(xs[ys==0,i],100,range=(0,15))
            counts1, bins1 = pl.histogram(xs[ys==1,i],100,range=(0,15))
        elif (i==7):
            counts0, bins0 = pl.histogram(xs[ys==0,i],100,range=(-1.5,1.))
            counts1, bins1 = pl.histogram(xs[ys==1,i],100,range=(-1.5,1.))	      
        else:
            counts0, bins0 = pl.histogram(xs[ys==0,i],100)
            counts1, bins1 = pl.histogram(xs[ys==1,i],100)
        pl.hold()
        pl.subplot(2,4,i+1)
        pl.plot(bins0[0:100],counts0,'r',bins1[0:100],counts1,'b')
        pl.title(feature_names[i])
    pl.tight_layout()
    pl.savefig("../out/{0}/{1}".format(WHICH_EXP,name),bbox_inches='tight')

def _temp_plot_(spk, ax, stim=0, yy=0):

	exid = np.where(spk[0] < ne)[0]
	inid = np.where(spk[0] >= ne)[0]

	htex = pl.histogram(spk[1][exid], bins=sim_time/10, range=(0, sim_time))
	htin = pl.histogram(spk[1][inid], bins=sim_time/10, range=(0, sim_time))

	hr = pl.histogram(spk[0], bins=n, range=(0, n)) 

	ax.plot(spk[1][exid]*dt, spk[0][exid], 'r.', markersize=mksz, label='Exc: '+ str(np.round(len(exid)/ne)) )
	ax.plot(spk[1][inid]*dt, spk[0][inid], 'b.', markersize=mksz, label='Inh: '+str(np.round(len(inid)/ne)) )

	ax.set_yticks([0, 99, 199, 299, ne-1, n-1])
	ax.set_yticklabels([])
	ax.set_ylim(0-10, n+10)
	ax.set_xlim([0-10, sim_time+10])

	ax.set_xticklabels([])	

	divider = make_axes_locatable(ax)	

	axHisty = divider.append_axes("right", size=.5, pad=0.1)
	#adjust_spines(axHisty,['left', 'bottom'], outward=0)

	axHisty.plot(hr[0]/(Ts), hr[1][0:-1], color='k', lw=2)
	
	if stim == 0: 
	   pl.text(.85, .5, str(np.round(len(exid)/ne / Ts,1)) +' Hz', transform = axHisty.transAxes, color='r')
	   pl.text(.85, .85, str(np.round(len(inid)/ni /Ts,1))+' Hz', transform = axHisty.transAxes, color='b')
	else:
	    pl.text(.9, .5, str(np.round(len(exid)/ne /Ts,1)) +' Hz', transform = axHisty.transAxes, color='r')
	    pl.text(.9, .85, str(np.round(len(inid)/ni /Ts,1))+' Hz', transform = axHisty.transAxes, color='b')

	axHisty.set_yticks([0, 99, 199, 299, ne-1, n-1])
	axHisty.set_yticklabels([])
	axHisty.set_xticks([0, 10])

	axHisty.set_ylim(0-10, n+10)

	axHistx = divider.append_axes("bottom", 1.2, pad=0.3)
	#adjust_spines(axHistx,['left', 'bottom'], outward=0)

	axHistx.plot(htex[1][0:-1], htex[0], color='r', lw=2, label='Exc')
	axHistx.plot(htin[1][0:-1], htin[0], color='b', lw=2, label='Inh')	
	
	axHistx.set_yticks([0, 50, 100, 150])
	axHistx.set_yticklabels([])
	
	if yy == 1:	
	   axHistx.set_xlabel('Time (ms)')
	   axHistx.set_ylabel('Population spike count')
	   axHistx.set_yticklabels([0, 50, 100, 150])
	   pl.legend(loc=1, frameon=False, prop={'size':12.5})
	   
	   axHisty.set_xlabel('Firing rate \n (spikes/s)', size=10)

def semiLogFracFound(i,o,**pltKwds):
    from pylab import histogram,semilogx,loglog
    f=i[(o>0) &(i>0)]
    l=i[(o<1) &(i>0)]

    d=i.sum()
    fNorm=f/d
    lNorm=l/d
    
    fHist=histogram(fNorm,logspace(-6,-2.0,30))
    lHist=histogram(lNorm,logspace(-6,-2.0,30))
    semilogx(fHist[1][1:],fHist[0].astype(float)/(fHist[0]+lHist[0]),**pltKwds)

def galaxy_stellar_massftn():
    gadget2msun=10.e10
    boxsize = 47.0
    max_mag=-16.
    min_mag = -23.
    nbins=14
    hubble_h = 0.7
    model2_folder = "/mnt/lustre/scratch/cs390/AHF_halos/cubepm_131212_6_1728_47Mpc_ext2/mergertrees/outputs/"
    nore_folder = "/mnt/lustre/scratch/cs390/AHF_halos/cubepm_131212_6_1728_47Mpc_ext2/mergertrees/outputs_nore/"
    snaplist_file = "/mnt/lustre/scratch/cs390/AHF_halos/cubepm_131212_6_1728_47Mpc_ext2/mergertrees/cubep3m_zlist_out"
    observe_folder="/mnt/lustre/scratch/cs390/codes/cubepm_131212_6_1728_47Mpc_ext2/observed_UVL/"
    firstfile = 0
    lastfile = 215

    filter = LGalaxyStruct.properties_used
    filter['DiskMass'] = True
    filter['BulgeMass'] = True

    file_prefix = "SA_z8.06"    
    (nTrees,nGals,nTreeGals,gal) = read_lgal.readsnap_lgal(model2_folder,file_prefix,firstfile,lastfile,filter)
    massf = gadget2msun*gal['DiskMass']+gadget2msun*gal['BulgeMass']
    mass = [i for i in massf if i > 10.e6]
    mass = numpy.log10(mass)
    stellarmass = pylab.histogram(mass)
    print stellarmass
    massftn_y = stellarmass[0]
    massftn_x = []
    for i in range(len(stellarmass[0])):
        massftn_x.append((stellarmass[1][i]+stellarmass[1][i+1])/2.)

    pylab.rc('text', usetex=True)
    fig = pylab.figure()
    ax = fig.add_subplot(111)
    ax.plot(massftn_x,massftn_y,'r-')


    file_prefix = "SA_z8.06"    
    (nTrees,nGals,nTreeGals,gal) = read_lgal.readsnap_lgal(nore_folder,file_prefix,firstfile,lastfile,filter)
    massf = gadget2msun*gal['DiskMass']+gadget2msun*gal['BulgeMass']
    mass = [i for i in massf if i > 10.e6]
    mass = numpy.log10(mass)
    stellarmass = pylab.histogram(mass)
    print stellarmass
    massftn_y = stellarmass[0]
    massftn_x = []
    for i in range(len(stellarmass[0])):
        massftn_x.append((stellarmass[1][i]+stellarmass[1][i+1])/2.)

    ax.plot(massftn_x,massftn_y,'b-')
    ax.set_yscale("log")
    pylab.show()

def plot_hists(nus=[143,353],
               map1_name=None,
               map2_name=None,
               maskname='wmap_temperature_kq85_analysis_mask_r10_9yr_v5.fits',
               nside=2048,
               fwhm=0.0,
              bins=100,normed=True,
              atol=1e-6, ymin=0.01, ymax=None,
              xmin=-0.001, xmax=0.005):

    if map1_name is None:
        map1_name = 'HFI_SkyMap_{}_2048_R2.02_full.fits'.format(nus[0])
    label1 = '{} GHz'.format(nus[0])
    if map2_name is None:
        map2_name = 'HFI_SkyMap_{}_2048_R2.02_full.fits'.format(nus[1])
    label2 = '{} GHz'.format(nus[1])
   
    map1 = prepare_map( map1_name, field=0,
                        maskname=maskname,
                        nside_out=nside, fwhm=fwhm )
    map2 = prepare_map( map2_name, field=0,
                        maskname=maskname,
                        nside_out=nside, fwhm=fwhm )

    y1,x1 = pl.histogram(map1[np.where(np.negative(np.isclose(map1,0.,atol=atol)))],
                       bins=bins,normed=normed)
    bin1 = (x1[:-1] + x1[1:]) / 2.

    y2,x2 = pl.histogram(map2[np.where(np.negative(np.isclose(map2,0.,atol=atol)))],
                       bins=bins,normed=normed)
    bin2 = (x2[:-1] + x2[1:]) / 2.
    #return bin1,y1,bin2,y2
        

    fig = plt.figure()
    ax = plt.gca()
    
    ax.semilogy(bin1, y1, lw=3, label=label1,color='red')
    ax.semilogy(bin2, y2, lw=3, label=label2,color='gray')
    ax.set_xlim(xmin=xmin,xmax=xmax)
    ax.set_ylim(ymin=ymin, ymax=ymax)

    #ax.set_yscale('log')
    
    ax.set_xlabel('$\mu K$', fontsize=20)
    ax.set_yticks([])
    
    plt.draw()
    plt.legend(frameon=False, fontsize=20)

    plt.savefig('pdfs_{}GHz_{}GHz_fwhm{:.3}rad.pdf'.format(nus[0],nus[1],fwhm))

def semiLogHistLostFound(i,o):
    from pylab import histogram,semilogx

    f=i[(o>0) &(i>0)]
    l=i[(o<1) &(i>0)]

    d=i.sum()
    fNorm=f/d
    lNorm=l/d
    
    fHist=histogram(fNorm,logspace(-6,-2.0,30))
    lHist=histogram(lNorm,logspace(-6,-2.0,30))
    
    semilogx(fHist[1][1:],fHist[0],label='out>0')
    semilogx(lHist[1][1:],lHist[0],label='out=0')

def spike_make_diagram(ts, gids, name):
    pylab.figure()
    color_marker = "."
    color_bar = "blue"
    color_edge = "black"
    ylabel = "Neuron ID"

    hist_binwidth = 5.0

    ax1 = pylab.axes([0.1, 0.3, 0.85, 0.6])
    pylab.plot(ts, gids, color_marker)
    pylab.ylabel(ylabel)
    pylab.xticks([])
    xlim = pylab.xlim()

    pylab.axes([0.1, 0.1, 0.85, 0.17])
    t_bins = numpy.arange(numpy.amin(ts), numpy.amax(ts), hist_binwidth)
    n, bins = pylab.histogram(ts, bins=t_bins)
    t_bins = t_bins[:-1]                        # FixMe it must work without cutting the end value
    num_neurons = len(numpy.unique(gids))
    heights = (1000 * n / (hist_binwidth * num_neurons))
    pylab.bar(t_bins, heights, width=hist_binwidth, color=color_bar, edgecolor=color_edge)
    pylab.yticks([int(a) for a in numpy.linspace(0.0, int(max(heights) * 1.1) + 5, 4)])
    pylab.ylabel("Rate (Hz)")
    pylab.xlabel("Time (ms)")
    pylab.xlim(xlim)
    pylab.axes(ax1)

    pylab.title('Spike activity')
    pylab.draw()
    pylab.savefig(path + name + ".png", dpi=dpi_n, format='png')
    pylab.close()

def logLogRatioFoundLost(i,o,**pltKwds):
    
    from pylab import histogram,semilogx,loglog

    f=i[(o>0) &(i>0)]
    l=i[(o<1) &(i>0)]

    d=i.sum()
    fNorm=f/d
    lNorm=l/d
    
    fHist=histogram(fNorm,logspace(-6,-2.0,30))
    lHist=histogram(lNorm,logspace(-6,-2.0,30))
    
    #semilogx(fHist[1][1:],fHist[0].astype(float)/lHist[0])
    loglog(fHist[1][1:],fHist[0].astype(float)/lHist[0],**pltKwds)

def dcf(X,Y,T,num_bin,noise_std):

	'''
	This function implements the discrete correlation function 
	(DCF) delay estimation method described in Edelson RA, Krolik JH (1988) "The discrete correlation function - A new method for 		analyzing unevenly sampled variability data." The Astrophysical Journal 333: 646-659.

	'''

	#obtain the delta Ts
	deltaT=T[:,None]-T[None,:]
	#iu1 = np.triu_indices(len(T),1)
	iu1 = np.triu_indices(len(T))
	hist, bin_edges=pb.histogram(np.abs(deltaT[iu1]),num_bin)
	cent=bin_edges[0:len(bin_edges)-1]+np.diff(bin_edges)*.5
	dcf=np.zeros(len(cent))
	sigx=np.var(X)
	sigy=np.var(Y)
	muX=np.mean(X)
	muY=np.mean(Y)
	for i in range(0,len(cent)):
		for j in range(0,len(T)):
			for k in range(j,len(T)):
				if i<len(cent)-1:
					if (np.abs(deltaT[j,k])>=bin_edges[i])&(np.abs(deltaT[j,k])<bin_edges[i+1]):
						dcf[i]+=((X[j]-muX)*(Y[k]-muY))/np.sqrt((sigx-noise_std**2)*(sigy-noise_std**2))
				elif i==len(cent)-1:
					if (np.abs(deltaT[j,k])>=bin_edges[i])&(np.abs(deltaT[j,k])<=bin_edges[i+1]):
						dcf[i]+=((X[j]-muX)*(Y[k]-muY))/np.sqrt((sigx-noise_std**2)*(sigy-noise_std**2))


	dcf[hist>0]=dcf[hist>0]/hist[hist>0]
	return cent[np.argmax(dcf)],dcf,cent

def testCollisionsE8(n,d=8):
    M = pylab.eye(8,8)

    S = [0.0]*n
    C = [0]*n
    #generate distances and buckets
    for i in range(n):
        p = [random() for j in xrange(d)]
        q = [p[j] + (gauss(0,1)/(d**.5)) for j in xrange(d)]
        S[i]=distance(p,q,d)
        C[i]= int(decodeE8(dot(p,M)) == decodeE8(dot(q,M)))
    
    ranges = pylab.histogram(S,30)[1]   
    bucketsCol = [0]*len(ranges)
    bucketsDis = [0]*len(ranges)

    #fill buckets with counts 
    for i in xrange(n):
        k = len(ranges)-1
        while S[i] < ranges[k]:k=k-1
        if C[i]:bucketsCol[k]=bucketsCol[k]+1
        else:bucketsDis[k] = bucketsDis[k]+1
    print bucketsDis
    print ranges
    pylab.plot(ranges,[float(bucketsCol[i])/(float(bucketsDis[i]+.000000000001))  for i in range(len(ranges))],color='purple') 

    def get_bcc_pz(self,filename_lenscat):

        if self.prob_z == None:


            # filename_lenscat = os.environ['HOME'] + '/data/BCC/bcc_a1.0b/aardvark_v1.0/lenscats/s2n10cats/aardvarkv1.0_des_lenscat_s2n10.351.fit'
            # filename_lenscat = os.environ['HOME'] + '/data/BCC/bcc_a1.0b/aardvark_v1.0/lenscats/s2n10cats/aardvarkv1.0_des_lenscat_s2n10.351.fit'
            lenscat = tabletools.loadTable(filename_lenscat)

            if 'z' in lenscat.dtype.names:
                self.prob_z , _  = pl.histogram(lenscat['z'],bins=self.grid_z_edges,normed=True)
            elif 'z-phot' in lenscat.dtype.names:
                self.prob_z , _  = pl.histogram(lenscat['z-phot'],bins=self.grid_z_edges,normed=True)

            if 'e1' in lenscat.dtype.names:
                self.sigma_ell = np.std(lenscat['e1'],ddof=1)

def plot_phases(in_file, plot_type, plot_log):
    flags = ['histogram','phases']
    plot_flag = 0
    log_flag = 0

    def no_log(x):
        return x

    fig = pylab.figure(1)
    ax = fig.add_subplot(111)

    try:
        img = spimage.sp_image_read(in_file,0)
    except:
        raise IOError("Can't read %s." % in_file)

    values = img.image.reshape(pylab.size(img.image))

    if plot_log:
        log_function = pylab.log
    else:
        log_function = no_log

    if plot_type == PHASES:
        hist = pylab.histogram(pylab.angle(values),bins=500)
        ax.plot((hist[1][:-1]+hist[1][1:])/2.0,log_function(hist[0]))
    elif plot_flag == HISTOGRAM:
        hist = pylab.histogram2d(pylab.real(values),pylab.imag(values),bins=500)
        ax.imshow(log_function(hist[0]),extent=(hist[2][0],hist[2][-1],-hist[1][-1],-hist[1][0]),interpolation='nearest')
    else:
        ax.plot(pylab.real(values),pylab.imag(values),'.')
    return fig

def old_spike_psth(data, t1_ms = -250., t2_ms = 0., bin_ms = 10):
  """Uses data format returned by get_spikes"""
  spike_time_ms = data['spike times ms']
  N_trials = data['trials']
  t2_ms = pylab.ceil((t2_ms - t1_ms) / bin_ms)*bin_ms + t1_ms
  N_bins = (t2_ms - t1_ms) / bin_ms
  
  if N_trials > 0:
    all_spikes_ms = pylab.array([],dtype=float)
    for trial in range(len(spike_time_ms)):
      if spike_time_ms[trial] is None:
        continue
      idx = pylab.find((spike_time_ms[trial] >= t1_ms) & 
                       (spike_time_ms[trial] <= t2_ms))
      all_spikes_ms = \
        pylab.concatenate((all_spikes_ms, spike_time_ms[trial][idx]))
    spike_n_bin, bin_edges = \
      pylab.histogram(all_spikes_ms, bins = N_bins, 
                      range = (t1_ms, t2_ms), new = True)

    spikes_per_trial_in_bin = spike_n_bin/float(N_trials) 
    spike_rate = 1000*spikes_per_trial_in_bin/bin_ms
  else:
    spike_rate = pylab.nan
  
  bin_center_ms = (bin_edges[1:] + bin_edges[:-1])/2.0

  return spike_rate, bin_center_ms

    def get_bcc_pz(self):

        if self.prob_z == None:


            filename_lenscat = os.environ['HOME'] + '/data/BCC/bcc_a1.0b/aardvark_v1.0/lenscats/s2n10cats/aardvarkv1.0_des_lenscat_s2n10.351.fit'
            lenscat = tabletools.loadTable(filename_lenscat)
            self.prob_z , _  = pl.histogram(lenscat['z'],bins=self.grid_z_edges,normed=True)

def calculate_activity_histogram(spikes, total_neurons, bin=0.1):
    """Calculates histogram and bins specifically for neurons.

    Bins are provided in seconds instead of milliseconds."""
    hist, bin_edges = pylab.histogram(
        [spike[0] for spike in spikes],
        bins=pylab.ceil(max([spike[0] for spike in spikes])/bin))
    bin_edges = pylab.delete(bin_edges, len(bin_edges)-1) / 1000.
    return [[float(i)/total_neurons for i in hist], bin_edges]

	def getTimeHistogramm(self, color, name):
		diffs = plt.diff(self.times)
		## compute standard histogram
		y,x = plt.histogram(diffs, bins=plt.linspace(diffs.min(), diffs.max(), 500))

		## notice that len(x) == len(y)+1
		## We are required to use stepMode=True so that PlotCurveItem will interpret this data correctly.
		curve = pg.PlotCurveItem(x, y, stepMode=True, fillLevel=0, pen=color, brush=color, name=name)
		return curve

def histeq(im, nbr_bins = 256):
	""" Histogram equalization of a grayscale image. """
	# get image histogram
	imhist, bins = pl.histogram(im.flatten(), nbr_bins, normed = True)
	cdf = imhist.cumsum()  # cumulative distribution function
	cdf = 255 * cdf / cdf[-1]  # normalize
	# use linear interpolation of cdf to find new pixel values
	im2 = pl.interp(im.flatten(), bins[:-1], cdf)
	return im2.reshape(im.shape)