def analyzeSound(self):
        """ highlights N first peaks in frequency diagram
        """
        # on recharge les données 
        data = self.data
        sample_freq = self.sample_freq
        from scipy.fftpack import fftfreq
        freq_vect = fftfreq(data.size) * sample_freq
        
        # on trouve les maxima
        y0 = abs(fft(data))
#        y1 = abs(fft(data[:, 1]))
        maxi0 = ((diff(sign(diff(y0))) < 0) & (y0[1:-1] > y0.max()/10.)).nonzero()[0] + 1 # local max
        # maxi1 = ((diff(sign(diff(y1))) < 0) & (y1[1:-1] > y1.max()/10.)).nonzero()[0] + 1 # local max
        
        # fréquence
        ax = self.main_figure.figure.add_subplot(212)
        ax.plot(freq_vect[maxi0], y0[maxi0], "o")
        # ax.plot(freq_vect[maxi1], y1[maxi1], "o")
        
        # annotations au dessus d'une fréquence de coupure
        fc = 100
        for point in maxi0[(freq_vect[maxi0] > fc).nonzero()][:self.ui.spinBox.value()]:
            plt.annotate("%.2f" % freq_vect[point], (freq_vect[point], y0[point]))
#        for point in maxi1[(freq_vect[maxi0] > fc).nonzero()][:self.ui.spinBox.value()]:
#            plt.annotate("%.2f" % freq_vect[point], (freq_vect[point], y1[point]))
        
        self.ui.main_figure.canvas.draw()

        def smooth_gamma(gamma=flat_gamma, knots=knots, tau=smoothing**-2):
            # the following is to include a "noise floor" so that level value
            # zero prior does not exert undue influence on age pattern
            # smoothing
            gamma = gamma.clip(pl.log(pl.exp(gamma).mean()/10.), pl.inf)  # only include smoothing on values within 10x of mean

            return mc.normal_like(pl.sqrt(pl.sum(pl.diff(gamma)**2 / pl.diff(knots))), 0, tau)

def calcAUC(data, y0, lag, mgr, asym, time):
    """
    Calculate the area under the curve of the logistic function
    using its integrated formula
    [ A( [A-y0] log[ exp( [4m(l-t)/A]+2 )+1 ]) / 4m ] + At
    """

    # First check that max growth rate is not zero
    # If so, calculate using the data instead of the equation
    if mgr == 0:
        auc = calcAUCData(data, time)
    else:
        timeS = time[0]
        timeE = time[-1]
        t1 = asym - y0
        #try:
        t2_s = py.log(py.exp((4 * mgr * (lag - timeS) / asym) + 2) + 1)
        t2_e = py.log(py.exp((4 * mgr * (lag - timeE) / asym) + 2) + 1)
        #except RuntimeWarning as rw:
            # Exponent is too large, setting to 10^3
        #    newexp = 1000
        #    t2_s = py.log(newexp + 1)
        #    t2_e = py.log(newexp + 1)
        t3 = 4 * mgr
        t4_s = asym * timeS
        t4_e = asym * timeE

        start = (asym * (t1 * t2_s) / t3) + t4_s
        end = (asym * (t1 * t2_e) / t3) + t4_e
        auc = end - start

    if py.absolute(auc) == float('Inf'):
        x = py.diff(time)
        auc = py.sum(x * data[1:])
    return auc

 def _phase_map(self):
     self.dphi = (2*P.pi * self.nhop * P.arange(self.nfft/2+1)) / self.nfft
     A = P.diff(P.angle(self.STFT),1) # Complete Phase Map
     U = P.c_[P.angle(self.STFT[:,0]), A - P.atleast_2d(self.dphi).T ]
     U = U - P.np.round(U/(2*P.pi))*2*P.pi
     self.dPhi = U
     return U

 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 calc_f_matrix(self):
     '''Calculate F-matrix for step iCSD method'''
     el_len = self.coord_electrode.size
     h_val = abs(pl.diff(self.coord_electrode)[0])
     self.f_matrix = pl.zeros((el_len, el_len))
     for j in xrange(el_len):
         for i in xrange(el_len):
             if i != 0:
                 lower_int = self.coord_electrode[i] - \
                     (self.coord_electrode[i] - \
                      self.coord_electrode[i - 1]) / 2
             else:
                 lower_int = pl.array([0, self.coord_electrode[i] - \
                                       h_val/2]).max()
             if i != el_len-1:
                 upper_int = self.coord_electrode[i] + \
                     (self.coord_electrode[i + 1] - \
                      self.coord_electrode[i]) / 2
             else:
                 upper_int = self.coord_electrode[i] + h_val / 2
             
             self.f_matrix[j, i] = si.quad(self.f_cylinder, a=lower_int, \
                         b=upper_int, args=(self.coord_electrode[j]), \
                         epsabs=self.tol)[0] + \
                 (self.cond - self.cond_top) / (self.cond + self.cond_top) *\
                 si.quad(self.f_cylinder, a=lower_int, b=upper_int, \
                         args=(-self.coord_electrode[j]), \
                         epsabs=self.tol)[0]

 def set_pdf(self, x, p, Nrl = 1000):
   """Generate the lookup tables. 
   x is the value of the random variate
   pdf is its probability density
   cdf is the cumulative pdf
   inversecdf is the inverse look up table
   
   """
   
   self.x = x
   self.pdf = p/p.sum() #normalize it
   self.cdf = self.pdf.cumsum()
   self.inversecdfbins = Nrl
   self.Nrl = Nrl
   y = pylab.arange(Nrl)/float(Nrl)
   delta = 1.0/Nrl
   self.inversecdf = pylab.zeros(Nrl)    
   self.inversecdf[0] = self.x[0]
   cdf_idx = 0
   for n in xrange(1,self.inversecdfbins):
     while self.cdf[cdf_idx] < y[n] and cdf_idx < Nrl:
       cdf_idx += 1
     self.inversecdf[n] = self.x[cdf_idx-1] + (self.x[cdf_idx] - self.x[cdf_idx-1]) * (y[n] - self.cdf[cdf_idx-1])/(self.cdf[cdf_idx] - self.cdf[cdf_idx-1]) 
     if cdf_idx >= Nrl:
       break
   self.delta_inversecdf = pylab.concatenate((pylab.diff(self.inversecdf), [0]))

def beat_track(x, feature=LogFrequencySpectrum, **kwargs):
    """
    Scheirer beat tracker. Use output of comb filter bank on filterbank
                           sub-bands to estimate tempo, and comb filter state
                           to track beats.
    inputs:
       x        - the audio signal or filename to analyze
       feature  - the feature class to use [LogFrequencySpectrum]
       **kwargs - parameters to the feature extractor [nbpo=1, nhop=441]
    outputs:
       z      - per-tempo comb filter summed outputs
       tempos - the range of tempos in z
       D      - the differentiated half-wave rectified octave-band filterbank
                outputs
    """
    kwargs.setdefault('nhop', 441)
    kwargs.setdefault('nbpo', 1)
    F = feature(x, **kwargs)
    frame_rate = F.sample_rate / float(F.nhop)
    D = diff(F.X, axis=1)
    D[where(D < 0)] = 0
    tempos = range(40, 200, 4)
    z = zeros((len(tempos), D.shape[0]))
    for i, bpm in enumerate(tempos):  # loop over tempos to test
        t = int(round(frame_rate * 60. / bpm))  # num frames per beat
        alpha = 0.5**(2.0/t)
        b = [1 - alpha]
        a = zeros(t)
        a[0] = 1.0
        a[-1] = alpha
        z[i, :] = lfilter(b, a, D).sum(1)  # filter and sum sub-band onsets
    return z, tempos, D

def remove_discontinuity(value, xgap=10, ygap=200):
    """
    Remove discontinuity (sudden jump) in a series of values.
    Written by Denis, developed for LLC Fringe Counts data.
    value : list or numpy.array
    xgap  : "width" of index of the list/array to adjust steps
    ygap  : threshold value to detect discontinuity
    """
    difflist = pl.diff(value)
    discont_index = pl.find(abs(difflist) > ygap)

    if len(discont_index) == 0:
        return value
    else:
        discont_index = pl.append(discont_index, len(difflist))

    # find indice at discontinuities
    discont = {"start": [], "end": []}
    qstart = discont_index[0]
    for i in range(len(discont_index) - 1):
        if discont_index[i + 1] - discont_index[i] > xgap:
            qend = discont_index[i]
            discont["start"].append(qstart - xgap)
            discont["end"].append(qend + xgap)
            qstart = discont_index[i + 1]

    # add offsets at discontinuities
    result = pl.array(value)
    for i in range(len(discont["end"])):
        result[0 : discont["start"][i]] += result[discont["end"][i]] - result[discont["start"][i]]

    # remove the median
    result = result - pl.median(result)
    return result

    def calc_k_matrix(self):
        '''Calculate the K-matrix used by to calculate E-matrices'''
        el_len = self.coord_electrode.size
        # expanding electrode grid
        z_js = pl.zeros(el_len+2)
        z_js[1:-1] = self.coord_electrode
        z_js[-1] = self.coord_electrode[-1] + \
            pl.diff(self.coord_electrode).mean()
        
        c_vec = 1./pl.diff(z_js)
        # Define transformation matrices
        c_jm1 = pl.matrix(pl.zeros((el_len+2, el_len+2)))
        c_j0 = pl.matrix(pl.zeros((el_len+2, el_len+2)))
        c_jall = pl.matrix(pl.zeros((el_len+2, el_len+2)))
        c_mat3 = pl.matrix(pl.zeros((el_len+1, el_len+1)))
        
        for i in xrange(el_len+1):
            for j in xrange(el_len+1):
                if i == j:
                    c_jm1[i+1, j+1] = c_vec[i]
                    c_j0[i, j] = c_jm1[i+1, j+1]
                    c_mat3[i, j] = c_vec[i]
        
        c_jm1[-1, -1] = 0
        
        c_jall = c_j0
        c_jall[0, 0] = 1
        c_jall[-1, -1] = 1
        
        c_j0 = 0
        
        tjp1 = pl.matrix(pl.zeros((el_len+2, el_len+2)))
        tjm1 = pl.matrix(pl.zeros((el_len+2, el_len+2)))
        tj0 = pl.matrix(pl.eye(el_len+2))
        tj0[0, 0] = 0
        tj0[-1, -1] = 0

        for i in xrange(1, el_len+2):
            for j in xrange(el_len+2):
                if i == j-1:
                    tjp1[i, j] = 1
                elif i == j+1:
                    tjm1[i, j] = 1
        
        # Defining K-matrix used to calculate e_mat1-3
        return (c_jm1*tjm1 + 2*c_jm1*tj0 + 2*c_jall + c_j0*tjp1)**-1 * 3 * \
            (c_jm1**2 * tj0 - c_jm1**2 * tjm1 + c_j0**2 * tjp1 - c_j0**2 * tj0)

def eye_sample_insert_interval(R):
  tt = R.data['Trials']['eyeXData']['Trial Time']  
  n_trials = len(tt)
  d_esii = pylab.array([],dtype=float)
  for tr in range(n_trials):
    d_esii = pylab.concatenate((d_esii,pylab.diff(tt[tr])))

  return d_esii

    def analyzeSound(self):
        """ highlights N first peaks in frequency diagram
        """
        # on recharge les données 
        data = self.data
        sample_freq = self.sample_freq
        from scipy.fftpack import fftfreq
        freq_vect = fftfreq(data.size) * sample_freq
        
        # on trouve les maxima
        y0 = abs(fft(data))
#        y1 = abs(fft(data[:, 1]))
        maxi0 = ((diff(sign(diff(y0))) < 0) & (y0[1:-1] > y0.max()/10.)).nonzero()[0] + 1 # local max
        # maxi1 = ((diff(sign(diff(y1))) < 0) & (y1[1:-1] > y1.max()/10.)).nonzero()[0] + 1 # local max
        
        # fréquence
#        ax = self.main_figure.figure.add_subplot(212)
#        ax.plot(freq_vect[maxi0], y0[maxi0], "o")
#        # ax.plot(freq_vect[maxi1], y1[maxi1], "o")
        
        # annotations au dessus d'une fréquence de coupure
        fc = 100.
        max_freq = []
        freq_vect_2 = freq_vect[maxi0]
        maxy0 = y0[maxi0]
        maxy0_list = maxy0.tolist()
        while np.size(max_freq) < 12 :
            F = np.argmax(maxy0_list)
            maxy0_list[F] = 0
            print(freq_vect_2[F])
            if np.abs(freq_vect_2[F]) > fc :
                max_freq.append(F)
            for i in range(0,4) :
                if (F+2-i) in range (1,len(maxy0_list)) and np.abs(freq_vect_2[F+2-i]-freq_vect_2[F]) < 4 :
                    maxy0_list[F+2-i] = 0
                    
        print(freq_vect_2[max_freq])
        for point in max_freq:
            plt.annotate("%.2f" % freq_vect_2[point], (freq_vect_2[point], maxy0[point]))
        ax = self.main_figure.figure.add_subplot(212)
        ax.plot(freq_vect_2[max_freq], maxy0[max_freq], "o")
        
#        for point in maxi1[(freq_vect[maxi0] > fc).nonzero()][:self.ui.spinBox.value()]:
#            plt.annotate("%.2f" % freq_vect[point], (freq_vect[point], y1[point]))
#        
        self.ui.main_figure.canvas.draw()

def testMonotonicIncrease(aRRay):
    if len(gr.find(gr.diff(aRRay)<0))>0:
        print 'The array is NOT monotonic.'
        x=0
    else: 
        print 'The array is non-decreasing.'
        x=1
    return x

 def unimodal_rate(f=rate, age_indices=age_indices, tau=1.e5):
     df = pl.diff(f[age_indices])
     sign_changes = pl.find((df[:-1] > NEARLY_ZERO) & (df[1:] < -NEARLY_ZERO))
     sign = pl.ones(len(age_indices)-2)
     if len(sign_changes) > 0:
         change_age = sign_changes[len(sign_changes)/2]
         sign[change_age:] = -1.
     return -tau*pl.dot(pl.absolute(df[:-1]), (sign * df[:-1] < 0))

def smooth(y,smoothBeta=smoothBeta):
    m = len(y)
    
    p = pl.diff(pl.eye(m),3).transpose()
    A = pl.eye(m)+smoothBeta*pl.dot(p.transpose(),p)
    
    smoothY = pl.solve(A, y)
    
    return smoothY

	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 mu_age_derivative_potential(mu_age=mu_age,
                                 increasing_a0=pl.clip(parameters['increasing']['age_start']-ages[0], 0, len(ages)),
                                 increasing_a1=pl.clip(parameters['increasing']['age_end']-ages[0], 0, len(ages)),
                                 decreasing_a0=pl.clip(parameters['decreasing']['age_start']-ages[0], 0, len(ages)),
                                 decreasing_a1=pl.clip(parameters['decreasing']['age_end']-ages[0], 0, len(ages))):
     mu_prime = pl.diff(mu_age)
     inc_violation = mu_prime[increasing_a0:increasing_a1].clip(-pl.inf, 0.).sum()
     dec_violation = mu_prime[decreasing_a0:decreasing_a1].clip(0., pl.inf).sum()
     return -1.e12 * (inc_violation**2 + dec_violation**2)

    def calc_e_matrices(self):
        '''Calculate the E-matrices used by cubic spline iCSD method'''
        el_len = self.coord_electrode.size
        ## expanding electrode grid
        z_js = pl.zeros(el_len+2)
        z_js[1:-1] = self.coord_electrode
        z_js[-1] = self.coord_electrode[-1] + \
            pl.diff(self.coord_electrode).mean()
        
        ## Define transformation matrices
        c_mat3 = pl.matrix(pl.zeros((el_len+1, el_len+1)))
        
        for i in xrange(el_len+1):
            for j in xrange(el_len+1):
                if i == j:
                    c_mat3[i, j] = 1./pl.diff(z_js)[i]

        # Get K-matrix
        k_matrix = self.calc_k_matrix()
        
        # Define matrixes for C to A transformation:
        # identity matrix except that it cuts off last element:
        tja = pl.matrix(pl.zeros((el_len+1, el_len+2)))
        # converting k_j to k_j+1 and cutting off last element:
        tjp1a = pl.matrix(pl.zeros((el_len+1, el_len+2))) 

        # C to A
        for i in xrange(el_len+1):
            for j in xrange(el_len+2):
                if i == j-1:
                    tjp1a[i, j] = 1
                elif i == j:
                    tja[i, j] = 1
        
        # Define spline coeffiscients
        e_mat0 = tja    
        e_mat1 = tja*k_matrix
        e_mat2 = 3 * c_mat3**2 * (tjp1a-tja) - c_mat3 * \
                (tjp1a + 2 * tja) * k_matrix
        e_mat3 = 2 * c_mat3**3 * (tja-tjp1a) + c_mat3**2 * \
                (tjp1a + tja) * k_matrix
        
        return e_mat0, e_mat1, e_mat2, e_mat3

	def calcFrequencies(self, times):
		self.times = times
		diffs = plt.diff(times)
		self.total_time = times[-1] - times[0]
		self.mean_period =  plt.mean(diffs) 
		self.max_period = diffs.max() 
		self.min_period = diffs.min()
		self.max_frequency = 1/ self.min_period
		self.min_frequency = 1 / self.max_period
		self.mean_frequency = 1 / self.mean_period

 def calc_f_inv_matrix(self):
     '''Calculate the inverse F-matrix for the standard CSD method'''
     h_val = abs(pl.diff(self.coord_electrode)[0])
     
     #Inner matrix elements  is just the discrete laplacian coefficients
     self.f_inv_matrix[0, 0] = -1
     for j in xrange(1, self.f_inv_matrix.shape[0]-1):
         self.f_inv_matrix[j, j-1:j+2] = [1., -2., 1.]
     self.f_inv_matrix[-1, -1] = -1
     
     self.f_inv_matrix = self.f_inv_matrix * -self.cond / h_val**2