def make_disp_map(mode, screen_res, fig_width, fig_height, fig_no):
         
      if mode == 'diagonal':
         max_x = int(pl.floor((screen_res[0]-fig_width)/20))
         max_y = int(pl.floor((screen_res[1]-fig_height)/20))
      elif mode == 'mosaic':
         max_x = int(pl.floor(screen_res[0]/fig_width))
         max_y = int(pl.floor(screen_res[1]/fig_height))
         
      disp_map = [] #[(0, 60, fig_width, fig_height)]

      top_offset = 60
      
      for y in range(0,max_y):
         if y == max_y-1:
            mx = max_x-1
         else:
            mx = max_x
         for x in range(0,mx):
            if mode=='mosaic':
               disp_map.append((x*fig_width, top_offset+y*fig_height, fig_width, fig_height))
            elif mode=='diagonal':
               disp_map.append((20*x*y, top_offset+15*y*x, fig_width, fig_height))
               
      return disp_map[np.mod(fig_no, max_x*max_y-1)]

def jetWoGn(reverse=False):
    """
    jetWoGn(reverse=False)
       - returning a colormap similar to cm.jet, but without green.
         if reverse=True, the map starts with red instead of blue.
    """
    m=18 # magic number, which works fine
    m0=pl.floor(m*0.0)
    m1=pl.floor(m*0.2)
    m2=pl.floor(m*0.2)
    m3=pl.floor(m/2)-m2-m1

    b_ = np.hstack( (0.4*np.arange(m1)/(m1-1.)+0.6, np.ones((m2+m3,)) ) )
    g_ = np.hstack( (np.zeros((m1,)),np.arange(m2)/(m2-1.),np.ones((m3,))) )
    r_ = np.hstack( (np.zeros((m1,)),np.zeros((m2,)),np.arange(m3)/(m3-1.)))

    r = np.hstack((r_,pl.flipud(b_)))
    g = np.hstack((g_,pl.flipud(g_)))
    b = np.hstack((b_,pl.flipud(r_)))

    if reverse:
        r = pl.flipud(r)
        g = pl.flipud(g)
        b = pl.flipud(b)

    ra = pl.linspace(0.0,1.0,m)

    cdict = {'red': zip(ra,r,r),
            'green': zip(ra,g,g),
            'blue': zip(ra,b,b)}

    return pl.matplotlib.colors.LinearSegmentedColormap('new_RdBl',cdict,256)

 def __init__(self, scanList, stateArray):
     if len(scanList)!=stateArray.shape[0]:
         raise Exception('number of scan and state should be the same')
     times = [scan.timestamp for scan in scanList]
     self.avgTime = times[int(pl.floor(len(times)/2))]
     #self.avgTime = pl.mean([scan.timestamp for scan in scanList])
     
     #transform the 3d coordinates of each scan
     #and numpy.vstack all the output m*3 array together
     
     #find average Lidar frame
     avgBodyState = stateArray[int(pl.floor(len(stateArray)/2))]
     #avgBodyState=pl.mean(stateArray, 0)
     
     
     w_R_avg_L, w_T_avg_L = self._bodyState2LidarState(avgBodyState)
     self.avgLidarState = self._matrix2State(w_R_avg_L, w_T_avg_L)
     
     transform = self._transformPointsFromBodyToAvgLidar                
     #from data points with transformation to avgState
     self.dataPoints = pl.vstack([transform(scan.dataArray, state, w_R_avg_L, w_T_avg_L)  for scan, state in zip(scanList, stateArray) if scan.hasValidData()])
     self.intensity = pl.vstack([scan.intensityArray for scan in scanList if scan.hasValidData()]).flatten()
     
     if self.dataPoints.shape[0]!=self.intensity.shape[0]:
         raise Exception('dist and intensity have different size')

 def _pvoc(self, X_hat, Phi_hat=None, R=None):
     """
     ::
       a phase vocoder - time-stretch
       inputs:
         X_hat - estimate of signal magnitude
         [Phi_hat] - estimate of signal phase
         [R] - resynthesis hop ratio
       output:
         updates self.X_hat with modified complex spectrum
     """
     N = self.nfft
     W = self.wfft
     H = self.nhop
     R = 1.0 if R is None else R
     dphi = (2*P.pi * H * P.arange(N/2+1)) / N
     print "Phase Vocoder Resynthesis...", N, W, H, R
     A = P.angle(self.STFT) if Phi_hat is None else Phi_hat
     phs = A[:,0]
     self.X_hat = []
     n_cols = X_hat.shape[1]
     t = 0
     while P.floor(t) < n_cols:
         tf = t - P.floor(t)            
         idx = P.arange(2)+int(P.floor(t))
         idx[1] = n_cols-1 if t >= n_cols-1 else idx[1]
         Xh = X_hat[:,idx]
         Xh = (1-tf)*Xh[:,0] + tf*Xh[:,1]
         self.X_hat.append(Xh*P.exp( 1j * phs))
         U = A[:,idx[1]] - A[:,idx[0]] - dphi
         U = U - P.np.round(U/(2*P.pi))*2*P.pi
         phs += (U + dphi)
         t += P.randn()*P.sqrt(PVOC_VAR*R) + R # 10% variance
     self.X_hat = P.np.array(self.X_hat).T

def dia1P_generateEnsamble(L,dilution,nSamp,distType,distPar,fNameBase):
	"""
	Generates nSamp samples of LxL lattice with given dilution. The 
	cluster sizes are generated according to the given distType:
	distType = 'P' : power law with exponent distPar
	distType = 'E' : exponential with exponent distPar
	
	The nodes are written to the file fNameBase.node
	The cluster size distribution is written to the file fNameBase.clusSize
	The cluster width distribution is written to the file fNameBase.clusWidth
	"""

	# Open files for writing
	nodeFile = open(fNameBase+'.node','w+')
	clusFile = open(fNameBase+'.clusSize','w+')
	widthFile = open(fNameBase+'.clusWidth','w+')

	# lists for storing the size and width distributions
	sizes = []
	widths = []

	nClus = 0		# Total number of clusters

	for sNum in range(nSamp):		# Working on the sNum^th instance
		# Generate cluster sizes
		cSizes = []
		if distType == 'P':	
			# Find the number of variables to draw
			nVar = pylab.floor(L*L*dilution*scipy.special.zeta(distPar,1)/scipy.special.zeta(distPar-1,1))
			cSizes = list(statUtils.discretePow(int(nVar),t=distPar).astype(int))
		elif distType == 'E':
			nVar = pylab.floor(L*L*dilution*(1 - pylab.exp(-distPar)))
			cSizes = list(statUtils.discreteExpn(int(nVar),l=distPar).astype(int))

		# Generate clusters
		clusters, n1, n2 = dia1P_generateDualClusters(L,cSizes)	
		# Write clusters to the file
		dia1P_writeClus(nodeFile,L,clusters)

		nClus = len(clusters)		# Update the total number of clusters

		# Update size and width distributions
		for clus in clusters:
			sz = dia1P_dualClusterSize(clus,L,statType='S')							# Calculate the size of the cluster
			wd = dia1P_dualClusterSize(clus,L,statType='W')							# Calculate the width of the cluster
			sizes.append(sz)
			widths.append(wd)
	
	# Write distributions to files
	for size in sizes:
		clusFile.write(str(size) + '\n')

	for width in widths:
		widthFile.write(str(width) +'\n')

	# Close files
	nodeFile.close()
	clusFile.close()
	widthFile.close()

def printTrack( fid , resized , frame , pos , sz ):
    p0 = [pos[0]-pylab.floor(sz[0]/2),pos[1]-pylab.ceil(sz[1]/2)]
    p1 = [pos[0]+pylab.floor(sz[0]/2),pos[1]+pylab.ceil(sz[1]/2)]

    if resized:
        p0 = [x*2 for x in p0]
        p1 = [x*2 for x in p1]

    fid.write(str(frame)+","+str(p0[1])+","+str(p0[0])+","+str(p1[1])+","+str(p1[0])+"\n")

def binvalues(bin_number,base):
    #returns binvalue, bin min, bin max for given bin number and base
    bin_max=float(base)
    bin_min=1.0
    for i in range(bin_number):
        bin_min*=base
        bin_max*=base
    bin_min=pl.floor(bin_min)
    bin_max=pl.floor(bin_max)
    if (bin_min!=bin_max):
        bin_min+=1.0
    return (pl.sqrt(bin_min*bin_max),bin_min,bin_max)

def wave2sram(waveA,waveB):#{{{
    r'''Construct sram sequence for waveform. This takes python array and converts it to sram data for hand off to the FPGA. Wave cannot be nddata - this is stupid...'''
    if not len(waveA)==len(waveB):
        raise Exception('Lengths of DAC A and DAC B waveforms must be equal.')
    dataA=[long(py.floor(0x1FFF*y)) for y in waveA] # Multiply wave by full scale of DAC
    dataB=[long(py.floor(0x1FFF*y)) for y in waveB]
    truncatedA=[y & 0x3FFF for y in dataA] # Chop off everything except lowest 14 bits
    truncatedB=[y & 0x3FFF for y in dataB]
    dacAData=truncatedA
    dacBData=[y<<14 for y in truncatedB] # Shift DAC B data by 14 bits, why is this done??
    sram=[dacAData[i]|dacBData[i] for i in range(len(dacAData))] # Combine DAC A and DAC B
    return sram#}}}

def choose_patches(IMAGES, L, batch_size=1000):
    sz = int(sqrt(L))
    imsz = shape(IMAGES)[0]
    num_images = shape(IMAGES)[2]
    BUFF = 4

    X = matrix(zeros([L,batch_size],'d'))
    for i in range(batch_size):
        j = int(floor(num_images * rand()))
        r = sz/2+BUFF+int(floor((imsz-sz-2*BUFF)*rand()))
        c = sz/2+BUFF+int(floor((imsz-sz-2*BUFF)*rand()))
        X[:,i] = reshape(IMAGES[r-sz/2:r+sz/2,c-sz/2:c+sz/2,j],[L,1])
    return X

def elapsed_to_abs_time(seconds, base_time=None):
    """
    Convert relative (elapsed) time in seconds to astronomical datetime.
    In:
        seconds : float, elapsed seconds
        base_time : datetime, recording start time (if available)
    Out:        
        datetime
    """
    if not base_time:
        base_time = dt.datetime(1900, 01, 01, 0, 0, 0,tzinfo=None)
    return base_time + dt.timedelta(seconds=pl.floor(seconds),\
        microseconds=(seconds - pl.floor(seconds))*(10**6))

def get_subwindow(im, pos, sz, cos_window):
    """
    Obtain sub-window from image, with replication-padding.
    Returns sub-window of image IM centered at POS ([y, x] coordinates),
    with size SZ ([height, width]). If any pixels are outside of the image,
    they will replicate the values at the borders.

    The subwindow is also normalized to range -0.5 .. 0.5, and the given
    cosine window COS_WINDOW is applied
    (though this part could be omitted to make the function more general).
    """

    if pylab.isscalar(sz):  # square sub-window
        sz = [sz, sz]

    ys = pylab.floor(pos[0]) \
        + pylab.arange(sz[0], dtype=int) - pylab.floor(sz[0]/2)
    xs = pylab.floor(pos[1]) \
        + pylab.arange(sz[1], dtype=int) - pylab.floor(sz[1]/2)

    ys = ys.astype(int)
    xs = xs.astype(int)

    # check for out-of-bounds coordinates,
    # and set them to the values at the borders
    ys[ys < 0] = 0
    ys[ys >= im.shape[0]] = im.shape[0] - 1

    xs[xs < 0] = 0
    xs[xs >= im.shape[1]] = im.shape[1] - 1
    #zs = range(im.shape[2])

    # extract image
    #out = im[pylab.ix_(ys, xs, zs)]
    out = im[pylab.ix_(ys, xs)]

    if debug:
        print("Out max/min value==", out.max(), "/", out.min())
        pylab.figure()
        pylab.imshow(out, cmap=pylab.cm.gray)
        pylab.title("cropped subwindow")

    #pre-process window --
    # normalize to range -0.5 .. 0.5
    # pixels are already in range 0 to 1
    out = out.astype(pylab.float64) / 255 - 0.5

    # apply cosine window
    out = pylab.multiply(cos_window, out)

    return out

    def as_sdms(self,decimals=0):
        min_val_size = len(str(int(floor(abs(self.lower_bound)*180/pi))))
        max_val_size = len(str(int(floor(abs(self.upper_bound)*180/pi))))
        deg_size=max(min_val_size, max_val_size)
        sign_char = '- +'[int(sign(self.value))+1]
        d_float   = abs(self.value)*180/pi
        d_int     = int(floor(d_float))
        m_float   = 60*(d_float - d_int)
        m_int     = int(floor(m_float))
        s_float   = 60*(m_float - m_int)
        s_int     = int(floor(s_float))
        
        frac_int    = int(floor(10**decimals*(s_float - s_int)+0.5))

        if frac_int >= 10**decimals:
            frac_int -= 10**decimals
            s_int +=1
        if s_int >= 60:
            s_int -= 60
            m_int += 1
        if m_int >= 60:
            m_int -= 60
            d_int += 1
        max_d = int(floor(self.upper_bound*180/pi+0.5))
        min_d = int(floor(self.lower_bound*180/pi+0.5))
        if d_int >= max_d and self.cyclical and not self.include_upper_bound:
            d_int -= (max_d-min_d)
        
        base_str = sign_char+str(d_int).rjust(deg_size,'0')+':'+str(m_int).rjust(2,'0')+':'+str(s_int).rjust(2,'0')
        if decimals is 0:
            return base_str
        else:
            return base_str+'.'+str(frac_int).rjust(decimals,'0')

def check(nelx,nely,rmin,x,dc):
    dcn=py.zeros((nely,nelx));
    for i in range(1,nelx+1):
      for j in range(1,nely+1):
          sumx=0.0 
          for k in range(py.maximum(i-py.floor(rmin),1),py.minimum(i+py.floor(rmin),nelx)+1):
              
              for l in range(py.maximum(j-py.floor(rmin),1),py.minimum(j+py.floor(rmin),nely)+1):
                fac = rmin-py.sqrt((i-k)**2+(j-l)**2)
                sumx = sumx+py.maximum(0,fac)
                dcn[j-1,i-1] = dcn[j-1,i-1] + py.maximum(0,fac)*x[l-1,k-1]*dc[l-1,k-1]
                 
          dcn[j-1,i-1] = dcn[j-1,i-1]/(x[j-1,i-1]*sumx)          
    
    return dcn   

    def _make_log_freq_map(self):
        """
        ::

            For the given ncoef (bands-per-octave) and nfft, calculate the center frequencies
            and bandwidths of linear and log-scaled frequency axes for a constant-Q transform.
        """
        fp = self.feature_params
        bpo = float(self.nbpo) # Bands per octave
        self._fftN = float(self.nfft)
        hi_edge = float( self.hi )
        lo_edge = float( self.lo )
        f_ratio = 2.0**( 1.0 / bpo ) # Constant-Q bandwidth
        self._cqtN = float( P.floor(P.log(hi_edge/lo_edge)/P.log(f_ratio)) )
        self._dctN = self._cqtN
        self._outN = float(self.nfft/2+1)
        if self._cqtN<1: print "warning: cqtN not positive definite"
        mxnorm = P.empty(self._cqtN) # Normalization coefficients        
        fftfrqs = self._fftfrqs #P.array([i * self.sample_rate / float(self._fftN) for i in P.arange(self._outN)])
        logfrqs=P.array([lo_edge * P.exp(P.log(2.0)*i/bpo) for i in P.arange(self._cqtN)])
        logfbws=P.array([max(logfrqs[i] * (f_ratio - 1.0), self.sample_rate / float(self._fftN)) 
                         for i in P.arange(self._cqtN)])
        #self._fftfrqs = fftfrqs
        self._logfrqs = logfrqs
        self._logfbws = logfbws
        self._make_cqt()

 def add_zeroline(current_data):
     from pylab import plot, legend, xticks, floor, xlim,ylim
     t = current_data.t
     #legend(('surface','topography'),loc='lower left')
     plot(t, 0*t, 'k')
     n = int(floor(t.max()/1800.)) + 2
     xticks([1800*i for i in range(n)],[str(0.5*i) for i in range(n)])

def plot_viz_of_stochs(vars, viz_func, figsize=(8,6)):
    """ Plot autocorrelation for all stochs in a dict or dict of dicts
    
    :Parameters:
      - `vars` : dictionary
      - `viz_func` : visualazation function such as ``acorr``, ``show_trace``, or ``hist``
      - `figsize` : tuple, size of figure
    
    """
    pl.figure(figsize=figsize)

    cells, stochs = tally_stochs(vars)

    # for each stoch, make an autocorrelation plot for each dimension
    rows = pl.floor(pl.sqrt(cells))
    cols = pl.ceil(cells/rows)

    tile = 1
    for s in sorted(stochs, key=lambda s: s.__name__):
        trace = s.trace()
        if len(trace.shape) == 1:
            trace = trace.reshape((len(trace), 1))
        for d in range(len(pl.atleast_1d(s.value))):
            pl.subplot(rows, cols, tile)
            viz_func(pl.atleast_2d(trace)[:, d])
            pl.title('\n\n%s[%d]'%(s.__name__, d), va='top', ha='center', fontsize=8)
            tile += 1

def displayData(X):
    print "Visualizing"
    m, n = X.shape
    width = round(sqrt(n))
    height = width
    display_rows = int(floor(sqrt(m)))
    display_cols = int(ceil(m/display_rows))

    print "Cell width:", width
    print "Cell height:", height    
    print "Display rows:", display_rows
    print "Display columns:", display_cols
        
    display = zeros((display_rows*height,display_cols*width))

    # Iterate through the training sets, reshape each one and populate
    # the display matrix with the letter matrixes.    
    for xrow in range(0, m):
        rowindex = divide(xrow, display_cols)
        columnindex = remainder(xrow, display_cols)
        rowstart = int(rowindex*height)
        rowend = int((rowindex+1)*height)
        colstart = int(columnindex*width)
        colend = int((columnindex+1)*width)
        display[rowstart:rowend, colstart:colend] = X[xrow,:].reshape(height,width).transpose()
         
    imshow(display, cmap=get_cmap('binary'), interpolation='none')
    
    # Show plot without blocking
    draw()    

def render_network(A):
    [L, M] = shape(A)
    sz = int(sqrt(L))
    buf = 1
    A = asarray(A)

    if floor(sqrt(M)) ** 2 != M:
        m = int(sqrt(M / 2))
        n = M / m
    else:
        m = int(sqrt(M))
        n = m

    array = -ones([buf + m * (sz + buf), buf + n * (sz + buf)], "d")

    k = 0
    for i in range(m):
        for j in range(n):
            clim = max(abs(A[:, k]))
            x_offset = buf + i * (sz + buf)
            y_offset = buf + j * (sz + buf)
            array[x_offset : x_offset + sz, y_offset : y_offset + sz] = reshape(A[:, k], [sz, sz]) / clim

            k += 1
    return array

 def _pvoc2(self, X_hat, Phi_hat=None, R=None):
     """
     ::
       alternate (batch) implementation of phase vocoder - time-stretch
       inputs:
         X_hat - estimate of signal magnitude
         [Phi_hat] - estimate of signal phase
         [R] - resynthesis hop ratio
       output:
         updates self.X_hat with modified complex spectrum
     """
     N, W, H = self.nfft, self.wfft, self.nhop
     R = 1.0 if R is None else R
     dphi = P.atleast_2d((2*P.pi * H * P.arange(N/2+1)) / N).T
     print "Phase Vocoder Resynthesis...", N, W, H, R
     A = P.angle(self.STFT) if Phi_hat is None else Phi_hat
     U = P.diff(A,1) - dphi
     U = U - P.np.round(U/(2*P.pi))*2*P.pi
     t = P.arange(0,n_cols,R)
     tf = t - P.floor(t)
     phs = P.c_[A[:,0], U] 
     phs += U[:,idx[1]] + dphi # Problem, what is idx ?
     Xh = (1-tf)*Xh[:-1] + tf*Xh[1:]
     Xh *= P.exp( 1j * phs)
     self.X_hat = Xh

def get_filterbanks(nfilt=20,nfft=512,samplerate=16000,lowfreq=0,highfreq=None):
    """Compute a Mel-filterbank. The filters are stored in the rows, the columns correspond
    to fft bins. The filters are returned as an array of size nfilt * (nfft/2 + 1)

    :param nfilt: the number of filters in the filterbank, default 20.
    :param nfft: the FFT size. Default is 512.
    :param samplerate: the samplerate of the signal we are working with. Affects mel spacing.
    :param lowfreq: lowest band edge of mel filters, default 0 Hz
    :param highfreq: highest band edge of mel filters, default samplerate/2
    :returns: A numpy array of size nfilt * (nfft/2 + 1) containing filterbank. Each row holds 1 filter.
    """
    highfreq= highfreq or samplerate/2
    assert highfreq <= samplerate/2, "highfreq is greater than samplerate/2"
    
    # compute points evenly spaced in mels
    lowmel = hz2mel(lowfreq)
    highmel = hz2mel(highfreq)
    melpoints = pylab.linspace(lowmel,highmel,nfilt+2)
    # our points are in Hz, but we use fft bins, so we have to convert
    #  from Hz to fft bin number
    bin = pylab.floor((nfft+1)*mel2hz(melpoints)/samplerate)

    fbank = pylab.zeros([nfilt,nfft/2+1])
    for j in xrange(0,nfilt):
        for i in xrange(int(bin[j]),int(bin[j+1])):
            fbank[j,i] = (i - bin[j])/(bin[j+1]-bin[j])
        for i in xrange(int(bin[j+1]),int(bin[j+2])):
            fbank[j,i] = (bin[j+2]-i)/(bin[j+2]-bin[j+1])
    return fbank