def timebetweenpeaks(dataMatrix, Sensor,Min = False,order=10):
    timedifference = []
    if Min==False:
        for Sensors in Sensor:
            peaks = (argrelmax(dataMatrix[:,Sensors],order=order))
            distance = []
            for p in range(0,len(peaks[0])-1):
                distance.append((peaks[0][p+1]-peaks[0][p]))
            number=0
            for d in distance:
                number +=d
            timedifference.append((number/len(distance)))

    else:
        for Sensors in Sensor:
            peaks = (argrelmax(dataMatrix[:,Sensors],order=order))
            distance = []
            for p in range(0,len(peaks[0])-1):
                distance.append((peaks[0][p+1]-peaks[0][p]))
            number=0
            for d in distance:
                number +=d
            timedifference.append((number/len(distance)))

    return timedifference

    def detect_jumps(self):
        wle = int(self.regressor_t/self.pitch_t_hop)
        isel = self.mag > np.max(self.mag)*self.mag_threshold
        tsel = self.t[isel]
        fsel = self.f0[isel]
        self.isel = isel
        #pl.plot(np.flatnonzero(isel),20*np.log10(m[isel]))

        le = linreg2_err(tsel, fsel, wleft=wle, wright=wle, use_l=True)
        #ax[0].plot(tsel,fsel)
        #ax[1].plot(tsel,le)
        
        imax = argrelmax(le)[0]
        lemax = le[imax]
        idx = imax[lemax > self.t_threshold]
        ijup = idx
        #ax[0].plot(tsel[idx],fsel[idx],'o')
        #ax[1].plot(tsel[idx],le[idx],'o')
        
        imin = argrelmax(-le)[0]
        lemin = le[imin]
        idx = imin[lemin < -self.t_threshold]
        ijdn = idx
        #ax[0].plot(tsel[idx],fsel[idx],'o')
        #ax[1].plot(tsel[idx],le[idx],'o')

        self.down_jump_indices = np.asarray(ijdn)
        self.up_jump_indices = np.asarray(ijup)
        self.down_jump_times = tsel[ijdn]
        self.up_jump_times = tsel[ijup]

    def on_change2(pt):
        #maxi=sp.argrelmax(normList)[0]
        #pt=maxi[pt]
        
        fig=plt.figure(figsize=(20,10))
        
        gs = gridspec.GridSpec(2, 2)
        ax1 = plt.subplot(gs[:, 0])
        ax2 = plt.subplot(gs[0,1])
        ax3 = plt.subplot(gs[1,1])
        
        #ax=plt.subplot(1,2,1)
        ax1.plot(f,normList)
        
        
        
        ax1.plot(f[pt],normList[pt],'ko')
        #ax1.text(f[pt],normList[pt],str(f[pt])+ 'Hz')
        string='f={:.3f} Hz\nMode={:.0f}'.format(f[pt],pt)
        ax1.text(0.05, 0.95, string, transform=ax1.transAxes, fontsize=14,
            verticalalignment='top')
        
        ax1.set_xscale('log')
        ax1.set_yscale('log')
        
        #ax=plt.subplot(1,2,2)
        idxMode=myDMD_Uy.getIdxforFrq(f[pt])
        mode=myDMD_Uy.getMode(idxMode)
        ax2.imshow(np.real(mode),vmin=vmin,vmax=vmax,interpolation='nearest')
        
        uy=np.array(np.real(mode)[iRow,:])
        uy_imag=np.array(np.imag(mode)[iRow,:])
        ax3.plot(uy)
        ax3.plot(uy_imag,'r')
        maxi=sp.argrelmax(uy)[0]
        mini=sp.argrelmin(uy)[0]
        exti=np.sort(np.r_[maxi,mini])
        
        maxi_imag=sp.argrelmax(uy_imag)[0]
        mini_imag=sp.argrelmin(uy_imag)[0]
        exti_imag=np.sort(np.r_[maxi_imag,mini_imag])        
        
        print np.diff(exti)
        ax3.scatter(maxi,uy[maxi],marker=2)
        ax3.scatter(mini,uy[mini],marker=3)
        ax3.scatter(maxi_imag,uy_imag[maxi_imag],marker=2)
        ax3.scatter(mini_imag,uy_imag[mini_imag],marker=3)       

        ax3.set_xlim([0,np.real(mode).shape[1]])
        gamma=0
        print 'n=',L/(np.diff(maxi)*dx)+gamma
        print 'n=',L/(np.diff(mini)*dx)+gamma
        print 'n=',L/(np.diff(exti)*dx*2.0)+gamma
        
        print 'n=',L/(np.diff(maxi_imag)*dx)+gamma
        print 'n=',L/(np.diff(mini_imag)*dx)+gamma
        print 'n=',L/(np.diff(exti_imag)*dx*2.0)+gamma

def getWinSize(im, padFactor=4):
    h, w = im.shape
    im = np.pad(im, ((0, h * (padFactor - 1)), (0, w * (padFactor - 1))), mode='constant')

    spec = abs(np.fft.fft2(im))

    horizIs = sig.argrelmax(spec[0, :])[0]
    vertIs = sig.argrelmax(spec[:, 0])[0]
    maxHorizI = max(horizIs, key=lambda i: spec[0, i])
    maxVertI = max(vertIs, key=lambda i: spec[i, 0])
    return round(float(im.shape[1]) / maxHorizI), round(float(im.shape[0]) / maxVertI)

def locxymax(nda, order=1, mode='clip') :
    """For 2-d or 3-d numpy array finds mask of local maxima in x and y axes (diagonals are ignored) 
       using scipy.signal.argrelmax and return their product.

       @param nda - input ndarray
       @param order - range to search for local maxima along each dimension
       @param mode - parameter of scipy.signal.argrelmax of how to treat the boarder
    """
    shape = nda.shape
    size  = nda.size
    ndim = len(shape)

    if ndim< 2 or ndim>3 :
        msg = 'ERROR: locxymax nda shape %s should be 2-d or 3-d' % (shape)
        sys.exit(msg)

    ext_cols = argrelmax(nda, -1, order, mode)
    ext_rows = argrelmax(nda, -2, order, mode)
    
    indc = np.array(ext_cols, dtype=np.uint16)
    indr = np.array(ext_rows, dtype=np.uint16)

    msk_ext_cols = np.zeros(shape, dtype=np.uint16)
    msk_ext_rows = np.zeros(shape, dtype=np.uint16)

    if ndim == 2 :
        icr = indc[0,:] 
        icc = indc[1,:]
        irr = indr[0,:]
        irc = indr[1,:]

        msk_ext_cols[icr,icc] = 1
        msk_ext_rows[irr,irc] = 1

    elif ndim == 3 :
        ics = indc[0,:] 
        icr = indc[1,:] 
        icc = indc[2,:]
        irs = indr[0,:]
        irr = indr[1,:]
        irc = indr[2,:]

        msk_ext_cols[ics,icr,icc] = 1
        msk_ext_rows[irs,irr,irc] = 1

    #print 'nda.size:',   nda.size
    #print 'indc.shape:', indc.shape
    #print 'indr.shape:', indr.shape
    
    return msk_ext_rows * msk_ext_cols

def extr(x):
    """Extract the indices of the extrema and zero crossings.
    :param x: input signal
    :type x: array-like
    :return: indices of minima, maxima and zero crossings.
    :rtype: tuple
    """
    m = x.shape[0]

    x1 = x[:m - 1]
    x2 = x[1:m]
    indzer = find(x1 * x2 < 0)
    if np.any(x == 0):
        iz = find(x == 0)
        indz = []
        if np.any(np.diff(iz) == 1):
            zer = x == 0
            dz = np.diff([0, zer, 0])
            debz = find(dz == 1)
            finz = find(dz == -1) - 1
            indz = np.round((debz + finz) / 2)
        else:
            indz = iz
        indzer = np.sort(np.hstack([indzer, indz]))

    indmax = argrelmax(x)[0]
    indmin = argrelmin(x)[0]

    return indmin, indmax, indzer

def get_envelops(x, t=None):
    """ Find the upper and lower envelopes of the array `x`.
    """
    if t is None:
        t = np.arange(x.shape[0])
    maxima = argrelmax(x)[0]
    minima = argrelmin(x)[0]

    # consider the start and end to be extrema

    ext_maxima = np.zeros((maxima.shape[0] + 2,), dtype=int)
    ext_maxima[1:-1] = maxima
    ext_maxima[0] = 0
    ext_maxima[-1] = t.shape[0] - 1

    ext_minima = np.zeros((minima.shape[0] + 2,), dtype=int)
    ext_minima[1:-1] = minima
    ext_minima[0] = 0
    ext_minima[-1] = t.shape[0] - 1

    tck = interpolate.splrep(t[ext_maxima], x[ext_maxima])
    upper = interpolate.splev(t, tck)
    tck = interpolate.splrep(t[ext_minima], x[ext_minima])
    lower = interpolate.splev(t, tck)
    return upper, lower

def decode(file_name):
    border.rotate(file_name)
    image = Image.open("temp.png")
    q = border.find("temp.png")
    ind = sp.argmin(sp.sum(q, 1), 0)
    up_left = q[ind, 0] + 2
    up_top = q[ind, 1] + 2
    d_right = q[ind+1, 0] - 3
    d_bottom = q[ind-1, 1] - 3

    box = (up_left, up_top, d_right, d_bottom)
    region = image.crop(box)
    h_sum = sp.sum(region, 0)
    m = argrelmax(sp.correlate(h_sum, h_sum, 'same'))
    s = sp.average(sp.diff(m))
    m = int(round(d_right - up_left)/s)
    if m % 3 != 0:
        m += 3 - m % 3
    n = int(round(d_bottom - up_top)/s)
    if n % 4 != 0:
        n += 4 - n % 4
    s = int(round(s))+1

    region = region.resize((s*m, s*n), PIL.Image.ANTIALIAS)
    region.save("0.png")
    pix = region.load()
    matrix = mix.off(rec.matrix(pix, s, m, n))
    str2 = hamming.decode(array_to_str(matrix))

    return hamming.bin_to_str(str2)

def elliptical_orbit_to_events(t, w):
    """
    Convert an orbit to MIDI events using Cartesian coordinates and rules.

    Parameters
    ----------
    t : array_like
    w : array_like
    midi_pool : array_like

    """

    loop = gd.classify_orbit(w)

    # apocenters
    x,y,z = w.T[:3]
    r = np.sqrt(x**2 + y**2 + z**2)
    apo = np.array([argrelmax(rr)[0] for rr in r])

    # get periods
    periods = []
    for i in range(w.shape[1]):
        if np.any(loop[i] == 1):
            w2 = gd.align_circulation_with_z(w[:,i], loop[i])

            R = np.sqrt(w2[:,0]**2 + w2[:,1]**2)
            phi = np.arctan2(w2[:,1], w2[:,0]) % (2*np.pi)
            z = w2[:,2]

            # loop
            T1 = gd.peak_to_peak_period(t, R)
            T2 = gd.peak_to_peak_period(t, phi)
            T3 = gd.peak_to_peak_period(t, z)

        else:
            # box
            T1 = gd.peak_to_peak_period(t, w[:,i,0])
            T2 = gd.peak_to_peak_period(t, w[:,i,1])
            T3 = gd.peak_to_peak_period(t, w[:,i,2])

        periods.append([T1,T2,T3])

    freqs = (2*np.pi / np.array(periods)) * 10000.

    delays = []
    notes = []
    for j in range(w.shape[0]):
        _no = []
        for i in range(w.shape[1]):
            if j in apo[i]:
                _no.append(freqs[i].tolist())

        if len(_no) > 0:
            delays.append(t[j])
            notes.append(np.unique(_no).tolist())

    delays = np.array(delays)
    notes = np.array(notes)

    return delays, notes

    def find_current_peak_position(self, data):
        """
        Finds the current average peak position.

        Parameters
        ----------
            data : numpy.ndarray
            A phase-corrected datacube.

        Returns
        -------
            peak_position : int
            The argument of the highest local maximum.
        """
        self.info('Finding current peak position.')
        data = data.sum(axis=2)
        data = data.sum(axis=1)
        data = np.where(data < 0.75 * data.max(), 0, data)
        peaks = signal.argrelmax(data, axis=0, order=5)[0]
        self.info('Encountered {:d} peaks: '.format(len(peaks)))

        peaks_values = data[peaks]
        max_peaks_arg = np.argmax(peaks_values)
        peak = peaks[max_peaks_arg]

        return peak

    def apocenter(self, type=np.mean):
        """
        Estimate the apocenter(s) of the orbit. By default, this returns
        the mean apocenter. To get, e.g., the minimum apocenter,
        pass in ``type=np.min``. To get all apocenters, pass in
        ``type=None``.

        Parameters
        ----------
        type : func (optional)
            By default, this returns the mean apocenter. To return all
            apocenters, pass in ``None``. To get, e.g., the minimum
            or maximum apocenter, pass in ``np.min`` or ``np.max``.

        Returns
        -------
        apo : float, :class:`~numpy.ndarray`
            Either a single number or an array of apocenters.
        """
        r = self.r
        max_ix = argrelmax(r, mode='wrap')[0]
        max_ix = max_ix[(max_ix != 0) & (max_ix != (len(r)-1))]

        if type is not None:
            return type(r[max_ix])
        else:
            return r[max_ix]

def par_find_peaks_by_chan(info):
    """
    Parameters
    ----------
    p_spect_array: numpy.ndarray
        An array with dimensions frequencies x channels
    frequencies: numpy.ndarray
        An array of the frequencies used
    std_thresh: float
        Threshold in number of standard deviations above the corrected power spectra to be counted as a peak

    Returns
    -------
    peaks_all_chans: numpy.ndarray with type bool
        An array of booleans the same shape as p_spect_array, specifying if there is a peak at a given frequency
        and electrode
    """

    p_spect_array = info[0]
    frequencies = info[1]
    std_thresh = info[2]

    peaks_all_chans = np.zeros(p_spect_array.shape).astype(bool)
    for i, chan_data in enumerate(p_spect_array.T):
        x = sm.tools.tools.add_constant(np.log10(frequencies))
        model_res = sm.RLM(chan_data, x).fit()
        peak_inds = argrelmax(model_res.resid)
        peaks = np.zeros(x.shape[0], dtype=bool)
        peaks[peak_inds] = True
        above_thresh = model_res.resid > (np.std(model_res.resid) * std_thresh)
        peaks_all_chans[:,i] = peaks & above_thresh
    return peaks_all_chans

def find_peaks(x, threshold=None, order=1):
    """Finds local maxima of a function.

    Args:
        x: A data vector.
        threshold: Local maxima under this value will be discarded.
            If threshold is None, the function will return only the global
            maximum.
            Defaut value os None.
        order: Number of samples before and after a data point that have
            to be smaller than the data point to be considered a local maximum.
            If 'threshold' is None, this argument has no effect.
            Default: 1.
    Returns:
        out: A list of local maxima, numpy array type.
    """
    if threshold is not None:
        event_peaks = signal.argrelmax(x, order=int(order))[0]
        if event_peaks.size > 0:
            return event_peaks[x[event_peaks] > threshold]
        return event_peaks
    else:
        if x.size > 0:
            return np.array([np.argmax(x)])
        return np.array([])

def _get_psp_list(bins, neuron_model, di_param, timestep, simtime):
    '''
    Return the list of effective weights from a list of NEST connection
    weights.
    '''
    nest.ResetKernel()
    nest.SetKernelStatus({"resolution":timestep})
    # create neuron and recorder
    neuron = nest.Create(neuron_model, params=di_param)
    vm = nest.Create("voltmeter", params={"interval": timestep})
    nest.Connect(vm, neuron)
    # send the spikes
    times = [ timestep+n*simtime for n in range(len(bins)) ]
    sg = nest.Create("spike_generator", params={'spike_times':times,
                                                'spike_weights':bins})
    nest.Connect(sg, neuron)
    nest.Simulate((len(bins)+1)*simtime)
    # get the max and its time
    dvm = nest.GetStatus(vm)[0]
    da_voltage = dvm["events"]["V_m"]
    da_times = dvm["events"]["times"]
    da_max_psp = da_voltage[ argrelmax(da_voltage) ]
    da_min_psp = da_voltage[ argrelmin(da_voltage) ]
    da_max_psp -= da_min_psp
    if len(bins) != len(da_max_psp):
        raise InvalidArgument("simtime too short: all PSP maxima are not in \
range")
    else:
        plt.plot(da_times, da_voltage)
        plt.show()
        return da_max_psp

def get_peaks(freq, spectrum, threshold, outdir, metadata):
    """
    This gets the frequency of the peaks of a signal.
    """
    # Smooth the data by using repeated mean smoothing.
    radius = 2
    for i in range(3):
        smooth_spec, _ = mean_smooth(spectrum, radius)
        freq = freq[radius:-radius+1]
        spectrum = spectrum[radius:-radius+1]

    # Get the peaks from the smoothed spectum.
    prad = 4
    peak_index = sig.argrelmax(smooth_spec, order=prad)[0]

    # Remove "peaks" that are only noise fluctuations.
    peaks = []
    for i in peak_index:
        lower = max(i - prad,  0)
        upper = min(i + prad + 1, len(smooth_spec))
        segment = smooth_spec[lower:upper] - smooth_spec[i]
        if abs(np.min(segment)) > threshold:
            peaks.append(i)

    # Frequencies and the spectra.
    freq_peaks = np.array([freq[i] for i in peaks])
    spec_peaks = np.array([spectrum[i] for i in peaks])

    # Create a plot.
    plt.plot(freq, spectrum, 'k')
    plt.plot(freq, smooth_spec, 'r')
    plt.plot(freq_peaks, spec_peaks, 'bo')
    plot_tools(metadata, outdir, 'peak_find')

    return freq_peaks

    def test_amgauss(self):
        """Test if the gaussian amplitude modulator works correctly."""
        time_center = 63
        n_points = 128
        spread = 10
        signal = am.amgauss(n_points, time_center, spread)
        # parameters of the underlying gaussian function of the form
        # f(x) = a * exp( (-(x - b) **2) / (2 * (c ** 2)))
        a, b, c = 1, time_center, spread / np.sqrt(2 * pi)
        # Integral of a Gaussian is a * c * sqrt( 2 * pi)
        integral = a * c * np.sqrt(2 * pi)
        self.assertAlmostEqual(integral, signal.sum())

        # Other miscellaneous properties of a Gaussian
        maximum = argrelmax(signal)
        self.assertEqual(len(maximum), 1)
        self.assertEqual(maximum[0][0], time_center - 1)
        self.assertAlmostEqual(signal[time_center - 1], 1.0)

        self.assert_is_monotonic_increasing(signal[:(time_center - 1)])
        self.assert_is_monotonic_decreasing(signal[(time_center - 1):])

        infpl1 = np.floor(b - c).astype(int) - 1
        infpl2 = np.floor(b + c).astype(int)
        self.assert_is_convex(signal[:infpl1])
        self.assert_is_concave(signal[infpl1:infpl2])
        self.assert_is_convex(signal[infpl2:])

    def spline_max_growth_rate(self, s, droplow=False):

        ### N.B.: set parameter of -2.3 for dropping low OD values from analysis - i.e., OD 0.1###
        if droplow: data = np.where(self.log_data < -2.3, 'nan', self.log_data)
        else: data = self.log_data
        interpolator = interpolate.UnivariateSpline(self.elapsed_time, data, k=4, s=s)  #k can be 3-5
        der = interpolator.derivative()

        # Get the approximation of the derivative at all points
        der_approx = der(self.elapsed_time)

        # Get the maximum
        self.maximum_index = np.argmax(der_approx)
        self.growth_rate = der_approx[self.maximum_index]
        self.doubling_time = np.log(2)/self.growth_rate
        self.time_of_max_rate = self.elapsed_time[self.maximum_index]

        # Get estimates of lag time and saturation time from 2nd derivative
        der2 = der.derivative()
        der2_approx = der2(self.elapsed_time)
        try: self.lag_index = signal.argrelmax(der2_approx)[0][0]  # find first max
        except: self.lag_index = 0
        if self.lag_index > self.maximum_index: self.lag_index = 0
        self.lag_time = self.elapsed_time[self.lag_index]
        self.lag_OD = self.raw_data[self.lag_index]
        minima = signal.argrelmin(der2_approx)[0]  # find first min after maximum_index
        which_min = bisect.bisect(minima, self.maximum_index)
        try: self.sat_index = minima[which_min]
        except: self.sat_index = len(self.elapsed_time) - 1
        self.sat_time = self.elapsed_time[self.sat_index]
        self.sat_OD = self.raw_data[self.sat_index]

        self.spline = interpolator(self.elapsed_time)
        self.intercept = self.log_data[self.maximum_index] - (self.growth_rate*self.time_of_max_rate) # b = y - ax
        self.fit_y_values = [((self.growth_rate * x) + self.intercept) for x in self.elapsed_time]  # for plotting

def manual_detect(x, times, ker_len, clip, rate):

    kernel = np.ones((int(ker_len))) / ker_len
    x_smoothed = convolve(x, kernel)

    boundaries = argrelmax(x_smoothed)[0]
    boundaries = np.append(boundaries, len(x)-1)
    boundaries = np.insert(boundaries, 0, 0)
    boundaries = times[boundaries]

    # Optionaly clip all boundaries that are too close apart
    if clip > 0:
        y = [boundaries[0]]
        i = 0
        for j in range(1, len(boundaries)):
            if boundaries[j]-boundaries[i] >= clip:
                boundaries[i:j] = np.mean(boundaries[i:j])
                i = j
            j += 1

        for bound in boundaries:
            if bound != y[-1]:
                y.append(bound)
        boundaries = np.array(y)

    return boundaries

    def extractResponse(self):

        stimStarts = (self.stimStartInds) / self.downSamplingFactor
        stimEnds = (self.stimEndInds) / self.downSamplingFactor
        samplingRateDown = self.vibrationSignalDown.sampling_rate

        self.stimAmps = []
        self.stimFreqs = []
        self.responseVTraces = []
        self.stimTraces = []
        for (stD, endD, st, end) in zip(stimStarts, stimEnds, self.stimStartInds, self.stimEndInds):

            stimDown = self.vibrationSignalDown[stD:endD + 1]
            stimDownFFT = np.fft.rfft(stimDown, n=2048)
            self.stimFreqs.append(np.argmax(np.abs(stimDownFFT)) * samplingRateDown / 2 / len(stimDownFFT))

            stimAS = self.vibrationSignal[st:end + 1]
            stim = stimAS.magnitude
            allAmps = stim[np.concatenate((argrelmin(stim)[0], argrelmax(stim)[0]))]

            self.stimAmps.append(np.abs(allAmps).mean() * self.vibrationSignal.units)

            self.responseVTraces.append(self.voltageSignal[st:end + 1])

            self.stimTraces.append((stimAS - np.mean(stimAS)))

def peakFind(data,l,pos,order=1,_even=False):
    sets = re.sub('x',str(pos),'Report: xm (m)')
    print sets    
    iterr = 1

    while not _even:                
        iterr += 1
        mean = int(l[2]*iterr)
        maxlen = []
        phases = {}
        framepeak = {}
        for frame in data.iteritems():
            name = frame[0]
            frame = frame[1][sets]
            frameMean = mo.rolling_sum(frame**3,mean,center=True).values
            maxpeak = sig.argrelmax(frameMean,order=order)
            maxlen.append(len(maxpeak[0]))
            framepeak[name] = frame.index[maxpeak]
            phase = np.zeros([framepeak[name].shape[0]-1])
            for i in xrange(phase.shape[0]):
                phase[i] = framepeak[name][i+1]-framepeak[name][i]
            phases[name] = phase
        try:   
            phases = pd.DataFrame(phases)
            framePeaks = pd.DataFrame(framepeak)
            _even = True
        except ValueError:
            print maxlen
            print '----- All of the peaks have not been found ------'
            if iterr is 200:
                _even = True
                print ' ------ max iterations reached. Check data length ------'
                sets,phases,framePeaks = [],[],[]
    
    return [sets,phases,framePeaks]