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 generate_cord():
    """
    TODO: Pass the parameters as input arguments.

    Parameters:
      - Fs       : sampling frequency
      - F0       : frequency of the notes forming chord
      - gain     : gains of individual notes in the chord
      - duration : duration of the chord in second
      - alpha    : attenuation in KS algorithm
    """
    Fs = 48000

    # D2, D3, F3, G3, F4, A4, C5, G5
    F0 = 440 * pylab.array(
        [2**-(31.0/12), 2**-(19.0/12), 2**-(16.0/12), 2**(-14.0/12),
         2**-(4.0/12), 1.0, 2**(3.0/12), 2**(10.0/12)])
    gain = [1.2, 3.0, 1.0, 2.2, 1.0, 1.0, 1.0, 3.5]
    duration = 4.0
    alpha = 0.9785

    # Number of samples in the chord.
    nbsample_chord = Fs * duration

    # This is used to correct alpha later, so that all the notes
    # decay together (with the same decay rate).
    first_duration = pylab.ceil(nbsample_chord / pylab.round_(Fs/F0[0]))

    # Initialization.
    chord = pylab.zeros(nbsample_chord)

    for i, f in enumerate(F0):
        print("Working on %g / %g" % (i+1, len(F0)))
        # Get M and duration parameter.
        current_M = pylab.round_(Fs/f)
        current_duration = pylab.ceil(nbsample_chord / current_M)

        # Correct current alpha so that all the notes decay together
        # (with the same decay rate)
        current_alpha = alpha ** (first_duration / current_duration)

        # Let Paul's high D on the bass ring a bit longer.
        if i == 1:
            current_alpha = current_alpha ** 0.8

        # Generate input and output of KS algorithm.
        x = pylab.rand(current_M)
        y = ks(x, current_alpha, int(current_duration))
        y = y[:int(nbsample_chord)]
        
        # Construct the chord by adding the generated note (with the
        # appropriate gain).
        chord = chord + gain[i] * y
        
    return Fs, duration, chord

def plot_spike_histogram(spikes, bin=0.1, total_neurons=None):
    print 'Plotting activity histogram.'
    global figure_dir
    pylab.figure(get_free_figure_number())
    pylab.clf()
    pylab.hist([spike[0] for spike in spikes],
               bins=pylab.ceil(max([spike[0] for spike in spikes])/bin))
    if total_neurons is not None:
        pylab.ylim([0, total_neurons])
    pylab.xlim([0., pylab.ceil(max([spike[0] for spike in spikes]))])
    pylab.xlabel('Time, ms')
    pylab.ylabel('Amount of active neurons')
    pylab.title('Average network activity')
    pylab.savefig(os.path.join(figure_dir, 'activity.png'))

def getAvgGreenTime(intergreen1, intergreen2):
    doc = libxml2.parseFile('tls.out')

    lNS = doc.xpathEval("count(/tls-states/tlsstate[@phase='0'])")
    lWE = doc.xpathEval("count(/tls-states/tlsstate[@phase='2'])")

    lIG1 = doc.xpathEval("count(/tls-states/tlsstate[@phase='1'])")
    lIG2 = doc.xpathEval("count(/tls-states/tlsstate[@phase='3'])")

    doc.freeDoc()

    greenNS = lNS / ceil((lIG1 / intergreen1))
    greenWE = lWE / ceil((lIG2 / intergreen2))

    return greenWE, greenNS

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 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 panel():
  global Ras,Rms,Ie,Ies
  global tstart
  global pan,t0,t1

  #Pulse start time
  tstart = 0.1

  #Ras = [1.,500.,20.] #Ohm*cm
  #Rms = [200.,10000.,200.] #Ohm*cm^2
  
  Ras = [5.0, 500.0, 50.0] #Ohm*cm
  Rms = [200.0, 10000.0, 1000.0] #Ohm*cm^2
  
  Ie = 4. #nA
  Ies = pl.c_[pl.ones_like(ns.t)]
  Ies[:pl.ceil(tstart/ns.dt)] = 0.

  t0 = 0.101
  t1 = 0.6

  pan = simulationpanel.SimulationPanel()
  pan.Move((640,600))
  pan.setdict(globals())
  pan.addcommand('sim()')
  pan.addcommand('plottao()')
  pan.addvar('Ras')
  pan.addvar('Rms')
  pan.addvar('Ie')
  pan.addvar('t0')
  pan.addvar('t1')

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 decodefft(finf,data, dropheights = False):
    #output: decoded data with the number of heights reduced
    #two variables are added to the finfo class:
    #deco_num_hei, deco_hrange
    #data must be arranged: 
    #    (channels,heights,times) (C-style, profs change faster)
    #fft along the entire(n=None) acquired heights(axis=1), stores in data
    num_chan = data.shape[0]
    num_ipps = data.shape[2]
    num_codes = finf.subcode.shape[0]
    num_bauds = finf.subcode.shape[1]
    NSA = finf.num_hei + num_bauds - 1
    uppower = py.ceil(py.log2(NSA))
    extra = int(2**uppower - finf.num_hei)
    NSA = int(2**uppower)
    fft_code = py.fft(finf.subcode,n = NSA,axis=1).conj()
    data = py.fft(data,n=NSA,axis=1) #n= None: no cropped data or padded zeros
    for ch in range(num_chan):
        for ipp in range(num_ipps):
            code_i = ipp % num_codes
            data[ch,:,ipp] = data[ch,:,ipp] * fft_code[code_i,:]
    data=py.ifft(data,n=NSA,axis=1) #fft along the heightsm
    if dropheights:
        return data[:,:-extra-(num_bauds-1),:]
    else:
        return data[:,:-extra,:]

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 gauss1(s, func):
    '''Construct a 1-D Gaussian mask''' 
    # for sufficient result use ceil(6s) by ceil(6s) for a gaussian filter
    # read: http://en.wikipedia.org/wiki/Gaussian_blur for more explaination     
    s = float(s)
    r = int(ceil(3 * s))
    n = int(ceil(6 * s) + 1)
    
    # n * n zero matrix of floats
    gaussFilter = zeros((n), dtype=float)

    # fill gaussFilter[] with gaussian values on x
    for x in range(n):
        gaussFilter[x] = func(s, x - r)
    
    return gaussFilter

 def imshow(self,shape,only=None):
     from pylab import subplot,imshow,cm,title,sqrt,ceil
     if only is None:
         only=list(range(len(self.weights)))
 
     L=len(only)
     c=ceil(sqrt(L))
     r=ceil(L/c)
 
     for i,idx in enumerate(only):
         w=self.weights[idx]
         w=w.reshape(shape)
         subplot(r,c,i+1)
         imshow(w,cmap=cm.gray,interpolation='nearest')
         title('PC %d' % (idx))
         

    def show_blocking(self, plotfile=''):
        '''Print out the blocking data and show a graph of the behaviour of the standard error with block size.
        
If plotfile is given, then the graph is saved to the specifed file rather than being shown on screen.'''

        # print blocking output
        # header...
        print '%-11s' % ('# of blocks'),
        fmt = '%-14s %-15s %-18s '
        header = ('mean (X_%i)', 'std.err. (X_%i)', 'std.err.err. (X_%i)')
        for data in self.data:
            data_header = tuple(x % (data.data_col) for x in header)
            print fmt % data_header,
        for key in self.covariance:
            str = 'cov(X_%s,X_%s)' % tuple(key.split(','))
            print '%-14s' % (str),
        for key in self.combination_stats:
            fmt = ['mean (X_%s'+self.combination+'X_%s)', 'std.err. (X_%s'+self.combination+'X_%s)']
            strs = tuple([s % tuple(key.split(',')) for s in fmt])
            print '%-16s %-18s' % strs,
        print
        # data
        block_fmt = '%-11i'
        fmt = '%-#14.12g %-#12.8e  %-#18.8e  '
        for s in range(len(self.data[0].stats)):
            print block_fmt % (self.data[0].stats[s].block_size),
            for data in self.data:
                print fmt % (data.stats[s].mean, data.stats[s].se, data.stats[s].se_error),
            for cov in self.covariance.itervalues():
                print '%+-#14.5e' % (cov[s]),
            for comb in self.combination_stats.itervalues():
                print '%-#16.12f %-#18.12e' % (comb[s].mean, comb[s].se),
            print

        # plot standard error 
        if PYLAB:
            # one sub plot per data set.
            nplots = len(self.data)
            for (i, data) in enumerate(self.data):
                pylab.subplot(nplots, 1, i+1)
                blocks = [stat.block_size for stat in data.stats]
                se = [stat.se for stat in data.stats]
                se_error = [stat.se_error for stat in data.stats]
                pylab.semilogx(blocks, se, 'g-', basex=2, label=r'$\sigma(X_{%s})$' % (data.data_col))
                pylab.errorbar(blocks, se, yerr=se_error, fmt=None, ecolor='g')
                xmax = 2**pylab.ceil(pylab.log2(blocks[0]+1))
                pylab.xlim(xmax, 1)
                pylab.ylabel('Standard error')
                pylab.legend(loc=2)
                if i != nplots - 1:
                    # Don't label x axis points.
                    ax = pylab.gca()
                    ax.set_xticklabels([])
            pylab.xlabel('# of blocks')
            if plotfile:
                pylab.savefig(plotfile)
            else:
                pylab.draw()
                pylab.show()

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 plot(self,only=None):
        from pylab import plot,subplot,sqrt,ceil,title
    
        weights=self.weights_xh.T
        
        if only is None:
            only=list(range(len(weights)))
    
        L=len(only)
        c=ceil(sqrt(L))
        r=ceil(L/c)


        for i,idx in enumerate(only):
            w=weights[idx]
            subplot(r,c,i+1)
            plot(w,'-o')
            title('Filter %d' % (idx))

    def zoom(self, event):
        newxlim = numpy.array(self.sp_iq.get_xlim())
        curxlim = numpy.array(self.xlim)
        if(newxlim[0] != curxlim[0] or newxlim[1] != curxlim[1]):
            #xmin = max(0, int(ceil(self.sample_rate*(newxlim[0] - self.position))))
            #xmax = min(int(ceil(self.sample_rate*(newxlim[1] - self.position))), len(self.iq))
            xmin = max(0, int(ceil(self.sample_rate*(newxlim[0]))))
            xmax = min(int(ceil(self.sample_rate*(newxlim[1]))), len(self.iq))

            iq = numpy.array(self.iq[xmin : xmax])
            time = numpy.array(self.time[xmin : xmax])

            iq_psd, freq = self.dopsd(iq)

            self.draw_psd(freq, iq_psd)
            self.xlim = numpy.array(self.sp_iq.get_xlim())

            draw()

 def zeroPadd(self,fbins):
     #zero padd the underlying tdData such that the fbins afterwards are fbins
     spac=1/self._tdData.dt/fbins
     actlen=self._tdData.getLength()
     nozeros=py.ceil(spac-actlen)
     self._tdData.zeroPaddData(nozeros)
     #leave the old bnds in place
     bnds=[min(self.getfreqs()),max(self.getfreqs())]
     zpd=self._calculatefdData(self._tdData)
     self.setFDData(self.getcroppedData(zpd,bnds[0],bnds[1]))

def visualize_priors(priors):
    from matplotlib import pyplot as plt
    N = len(priors)
    m = round(sqrt(N*3/4.))
    n = ceil(N/m)
    fig = plt.figure()
    for i, var in enumerate(priors):
        ax = fig.add_subplot(n,m,i)
        ax = ax.hist(priors[var].rvs(1000), 20)
    plt.show()

    def zoom(self, event):
        newxlim = numpy.array(self.sp_iq.get_xlim())
        curxlim = numpy.array(self.xlim)
        if(newxlim[0] != curxlim[0] or newxlim[1] != curxlim[1]):
            self.xlim = newxlim
            #xmin = max(0, int(ceil(self.sample_rate*(self.xlim[0] - self.position))))
            #xmax = min(int(ceil(self.sample_rate*(self.xlim[1] - self.position))), len(self.iq))
            xmin = max(0, int(ceil(self.sample_rate*(self.xlim[0]))))
            xmax = min(int(ceil(self.sample_rate*(self.xlim[1]))), len(self.iq))

            iq = self.iq[xmin : xmax]
            time = self.time[xmin : xmax]

            iq_fft = self.dofft(iq)
            freq = self.calc_freq(time, self.sample_rate)

            self.plot_fft[0].set_data(freq, iq_fft)
            self.sp_fft.axis([freq.min(), freq.max(),
                              iq_fft.min()-10, iq_fft.max()+10])

            draw()

def show_images(images,which_images=None,max_images=None):
    from pylab import imshow,subplot,sqrt,ceil,title,cm,gca
    from random import shuffle
    
    if which_images is None:
        which_images=list(range(len(images.data)))
        
    if isinstance(which_images[0],str):  # target names
        which_names=which_images
        which_images=[]
        for idx in range(len(images.data)):
            name=images.target_names[images.targets[idx]]
            if name in which_names:
                which_images.append(idx)
        
    if not max_images is None:
        shuffle(which_images)
        which_images=which_images[:max_images]
    
        
        
    if not which_images:
        raise ValueError("No images selected")
        
    L=len(which_images)
    c=ceil(sqrt(L))
    r=ceil(L/c)

    for i,idx in enumerate(which_images):
        im=images.data[idx]
        name=images.target_names[images.targets[idx]]
        subplot(r,c,i+1)
        imshow(im,interpolation='nearest',cmap=cm.gray)
        title(name)
        
        if i<(L-c):
            gca().set_xticklabels([])
            
        if i%c!=0:
            gca().set_yticklabels([])