def WidthEstimate1D(inList, method="interpolate"):
    scales = np.zeros(len(inList))
    for idx, y in enumerate(inList):
        x = fft.fftfreq(len(y)) * len(y) / 2.0
        if method == "interpolate":
            minima = (argrelmin(y))[0]
            if minima[0] > 1:
                interpolator = interp1d(y[0 : minima[0]], x[0 : minima[0]])
                scales[idx] = interpolator(np.exp(-1))
        if method == "fit":
            g = models.Gaussian1D(amplitude=y[0], mean=[0], stddev=[10], fixed={"amplitude": True, "mean": True})
            fit_g = fitting.LevMarLSQFitter()
            minima = (argrelmin(y))[0]
            if minima[0] > 1:
                xtrans = (np.abs(x) ** 0.5)[0 : minima[0]]
                yfit = y[0 : minima[0]]
            else:
                xtrans = np.abs(x) ** 0.5
                yfit = y
            output = fit_g(g, xtrans, yfit)
            scales[idx] = np.abs(output.stddev.value[0]) * (2 ** 0.5)
    #             expmod = Model(Exponential1D)
    #             pars = expmod.make_params(amp=y[0],scale=5.0)
    #             pars['amp'].vary = False
    #             result = expmod.fit(y,x=x,params = pars)
    #             scales[idx] = result.params['scale'].value
    return scales

def cyl_orbit_to_events(t, w, midi_pool_hi, midi_pool_lo, time_resolution=None):
    """
    Convert an orbit to MIDI events using cylindrical coordinates and rules.

    For cylindrical orbits, crossing the disk midplane (x-y plane) triggers a
    high note. Crossing the x-z plane triggers a low note. The pitch of the note
    is set by the cylindrical radius at the time of either crossing. Smaller
    radius triggers a higher pitch note.

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

    """

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

    # set amplitudes from size of z oscillations
    all_amps = np.abs(z).max(axis=0) / 10.

    # variable length arrays
    phi_cross = np.array([argrelmin(pphi)[0] for pphi in phi.T])
    z_cross = np.array([argrelmin(zz**2)[0] for zz in z.T])

    # quantize R orbit
    RR = np.sqrt(R)
    nbins_hi = len(midi_pool_hi)
    q_R_hi = quantize(RR, nbins=nbins_hi, min=RR.max(), max=RR.min())
    nbins_lo = len(midi_pool_lo)
    q_R_lo = quantize(RR, nbins=nbins_lo, min=RR.max(), max=RR.min())

    delays = []
    notes = []
    amps = []
    for j in range(w.shape[0]):
        _no = []
        _amps = []
        for i in range(w.shape[1]):
            if j in z_cross[i]:
                _no.append(midi_pool_hi[q_R_hi[j,i]])
                _amps.append(all_amps[i])

            if j in phi_cross[i]:
                _no.append(midi_pool_lo[q_R_lo[j,i]])
                _amps.append(all_amps[i])

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

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

    return delays, notes, amps

    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 elliptical_orbit_to_events2(t, w, x_pool, y_pool, z_pool):
    """
    Convert an orbit to MIDI events using Cartesian coordinates and rules.

    For Cartesian orbits...

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

    """

    x,y,z = w.T[:3]

    # quantize the periods and map on to notes
    x_cross = np.array([argrelmin(xx**2)[0] for xx in x])
    y_cross = np.array([argrelmin(yy**2)[0] for yy in y])
    z_cross = np.array([argrelmin(zz**2)[0] for zz in z])

    r_x = np.sqrt(y**2 + z**2)
    r_y = np.sqrt(x**2 + z**2)
    r_z = np.sqrt(x**2 + y**2)

    q_r_x = quantize(np.sqrt(r_x), nbins=len(x_pool),
                     max=np.sqrt(r_x).min(), min=np.sqrt(r_x).max())
    q_r_y = quantize(np.sqrt(r_y), nbins=len(y_pool),
                     max=np.sqrt(r_y).min(), min=np.sqrt(r_y).max())
    q_r_z = quantize(np.sqrt(r_z), nbins=len(z_pool),
                     max=np.sqrt(r_z).min(), min=np.sqrt(r_z).max())

    delays = []
    notes = []
    for j in range(w.shape[0]):
        _no = []
        for i in range(w.shape[1]):
            if j in x_cross[i]:
                _no.append(x_pool[q_r_x[i,j]])

            if j in y_cross[i]:
                _no.append(y_pool[q_r_y[i,j]])

            if j in z_cross[i]:
                _no.append(z_pool[q_r_z[i,j]])

        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_peaks(self, xdata, ydata, domain, std_dev=11):
        from scipy.ndimage.filters import gaussian_filter
        from scipy.signal import argrelmin
        import numpy as np

        xpeaks, ypeaks = self.choose_domain(xdata, ydata, domain)

        ygauss = gaussian_filter(ypeaks, std_dev)

        x_peak_coord = xpeaks[argrelmin(ygauss)[0]]
        y_peak_coord = ypeaks[argrelmin(ygauss)[0]]

        return x_peak_coord, y_peak_coord

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 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 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 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 locate_crossings(s, spectrum):
    ''' '''

    N = spectrum.shape[1]   # number of traced states
    deltas = [spectrum[:,k+1]-spectrum[:,k] for k in range(N-1)]
    deltas = np.array(deltas)

    print(deltas.shape)
    inds, xs = argrelmin(deltas, axis=1, order=1)
    xings = sorted(zip(xs, inds), key=lambda x: (x[0], -x[1]))

    # causal filter
    n_active = 0
    causal_xings = []
    for x, n in xings:
        if n <= n_active:
            causal_xings.append((x, n))
            if n == n_active:
                n_active += 1

    print(causal_xings)

    for x, n in causal_xings:
        analyse_crossing(deltas[n,:], x)


    plt.figure('Deltas')
    plt.plot(s, deltas.T[:,:5], 'x')
    plt.show()

    def ground_check(self, occ, show=False):

        try:
            assert hasattr(self, 's') and hasattr(self, 'spectrum')
        except AssertionError:
            print('Spectrum has not yet been solved...')
            return None

        delta = self.spectrum[:,1]-self.spectrum[:,0]
        xs = argrelmin(delta)[0]

        # accept the minimum gap and all sufficients small gaps
        dmin = min(delta[x] for x in xs)

        if show:
            plt.plot(delta)
            plt.show()
        xs = [x for x in xs if delta[x] < GAP_MAX or delta[x]<1.1*dmin]
        if len(xs) != 1:
            return [None]

        # only one gap, extract prob
        dx = int(self.nsteps*.05)
        popt, perr = analyse_crossing(delta, xs[0], dx, show=show)

        print(popt)

        if popt is None:
            return [None,]

        prob = occ[0]*1./np.sum(occ)

        params = [(prob, popt[1], abs(popt[2])/self.nsteps)]
        return params

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 get_pls_min_index(self, pls_cutoff_frequency=PLS_SUPER_CUTOFF_FREQUENCY):
        # 3.5Hzでスムージングした曲線の最小値の+-200ms以内に最小値が存在する
        MIN_POINTS_WINDOW = 200 #+-200[ms]以内
        length = self.data.shape[0]
        num = length / FFT_WINDOW_NUM + 1
        pls_smooth = []
        logging.debug(num)
        for i in np.arange(num):
            logging.debug(i)
            if i == num - 1:
                pls_fftfreq = np.fft.fftfreq(length - FFT_WINDOW_NUM * (num - 1), 1.0 / SAMPLING_FREQUENCY)
                pls_fft = np.fft.fft( self.data[PLS_NORMALIZED].values[FFT_WINDOW_NUM * i:])
            else:
                pls_fftfreq = np.fft.fftfreq( FFT_WINDOW_NUM, 1.0/SAMPLING_FREQUENCY)
                pls_fft = np.fft.fft( self.data[PLS_NORMALIZED].values[FFT_WINDOW_NUM * i:FFT_WINDOW_NUM * (i + 1) ] )

            pls_fft[ (pls_fftfreq > pls_cutoff_frequency) | (pls_fftfreq < -pls_cutoff_frequency) ] = 0
            pls_smooth.extend( np.real(np.fft.ifft(pls_fft)) )
        self.data['PLS(smooth)'] = pls_smooth
        min_point_indexes_of_smooth = signal.argrelmin(self.data['PLS(smooth)'].values)[0]
        min_indexes = []
        for index in min_point_indexes_of_smooth:
            if index > 100:
                min_index = self.data[PLS_NORMALIZED].iloc[index - 100:index + 100].argmin()
            else:
                min_index = self.data[PLS_NORMALIZED].iloc[0:index + 100].argmin()
            min_indexes.append(min_index)
        min_indexes = pd.unique(min_indexes)
        self.data[PLS_MIN] = self.data[PLS_NORMALIZED][min_indexes]
        logging.debug(self.data[PLS_MIN])
        return self.data

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

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

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

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

def import_idl(filename = '../../../Data/K16/kic126corr_n.sav', num = 400):
	idlfile = rs(filename);
	ident = np.array(idlfile['cont']).astype('int');
	flux = np.array(idlfile['flux']).astype('float64');
	cadence = np.array(idlfile['time']).astype('float64');

	for a in np.unique(ident): #data quarter
		mean = np.average(flux[idlfile['cont'] == a]);
		data = [cadence[idlfile['cont'] == a], flux[idlfile['cont'] == a]];
		print len(data[0]), len(data[1]);
		k =  movavg_final(data, num);
		
		flux[np.where(idlfile['cont'] == a)[0]] -= k.astype('float64');

		#gp2(np.array([cadence[idlfile['cont'] == a], flux[idlfile['cont'] == a]]), block_size = 2000)[2];
		#flux[idlfile['cont'] == a] -= mean;
		flux[idlfile['cont'] == a] /= mean;

	arm = argrelmin(flux)[0];
	arm = arm[flux[arm] < -0.005]
	for u in arm:
		fluxbase = flux[max(0,u - int(num)):min(len(flux), u + int(num))];
		fluxbase_mean = np.average(fluxbase[fluxbase > 0])
		flux[max(0,u - int(num)):min(len(flux), u + int(num))] -= fluxbase_mean;

	return cadence, flux;

 def argers(gamma):
     inds=list(argrelmin(absolute(-gamma/f0*(sin(2*X)-2*X)/(2*X**2)-X/(Np*pi)))[0])
     data=[f0,]
     if len(inds)>0:
         data.extend(f[inds])
     else:
         data=[f0,f0,f0]
     return data

def compute_mins(xlocs, yvals, window_size = 10):
    yval_locs = argrelmin(yvals, order = window_size)[0]
    if len(yval_locs) == 0:
        return [];
    minPoints = xlocs[yval_locs]
    list1 = []
    for i in minPoints:
        list1.append(i)
    return list1

def getrelmaximas(dataMatrix, Sensor,Min = False):
    extrema = []
    if Min==False:
        for Sensors in Sensor:
            extrema.append(dataMatrix[argrelmax(dataMatrix[:,Sensors],order=10),Sensors][0])
    else:
        for Sensors in Sensor:
            extrema.append(dataMatrix[argrelmin(dataMatrix[:,Sensors],order=10),Sensors][0])
    return pd.DataFrame.transpose(pd.DataFrame(extrema))

 def test_residue(self):
     """Test the residue of the emd output."""
     signal = np.sum([self.trend, self.mode1, self.mode2], axis=0)
     decomposer = EMD(signal, t=self.ts)
     imfs = decomposer.decompose()
     n_imfs = imfs.shape[0]
     n_maxima = argrelmax(imfs[n_imfs - 1, :])[0].shape[0]
     n_minima = argrelmin(imfs[n_imfs - 1, :])[0].shape[0]
     self.assertTrue(max(n_maxima, n_minima) <= 2)

def min_data(x_data, y_data, width, no_peaks):
    min_ind = signal.argrelmin(y_data, np.array(width))
    #print peakind
    #plt.show()
    x_mins = map(lambda x: x_data[x], min_ind)
    #print time_peaks
    y_mins = map(lambda x: y_data[x], min_ind)
    #print p_peaks
    return x_mins, y_mins