def centroid(xarr, yarr, kern=default_kernal, mask=None, mode="same"):
    """Find the centroid of a line following a similar algorithm as
        the center1d algorithm in IRAF.   xarr and yarr should be an area
        around the desired feature to be centroided.  The default kernal
        is used if the user does not specific one.

        The algorithm solves for the solution to the equation

        ..math:: \int (I-I_0) f(x-x_0) dx = 0

        returns xc
    """
    if len(yarr) < len(kern):
        raise SpectrographError("Array has to be larger than kernal")

    if mask is not None:
        # catch the fact that at the edges it
        if mask.sum() < len(default_kernal):
            warr = np.convolve(yarr, kern, mode=mode)
            xc = np.interp(0, warr[mask], xarr[mask])
            return xc
        else:
            yarr = yarr[mask]
            xarr = xarr[mask]

    # convle the input array with the default kernal
    warr = np.convolve(yarr, kern, mode=mode)

    # interpolate the results
    xc = np.interp(0, warr, xarr)

    return xc

    def __init__(self, 
                 simulator,
                 n = 60 * 60 * 2,  # 1 hour
                 stdev_coef = 1.5,
                 days = 5,
                 trading_hours_per_day = 3.75,
                 tick_size = 1,
                 smooth_n = 3
                 ):
        self.__dict__.update(locals())
        del self.self
        
        self.price = self.simulator.history['last_price'].astype('int')
        self.vol = self.simulator.history['vol']
        self.m = int(days * trading_hours_per_day * 60 * 60 * 2)
        self.w = np.ones(self.n + 1) / float(self.n + 1)
        self.y = self.price - np.convolve(self.price, self.w, mode='same')
        denorm = np.sqrt(np.convolve(self.y**2, self.w, mode='same'))
        denorm_mean = denorm.mean()
        denorm2 = (denorm > denorm_mean) * denorm + (denorm <= denorm_mean) * denorm_mean
#        self.ny = self.y
#        self.ny = (denorm > 0) * self.y / denorm
        self.ny = (denorm2 > 0) * self.y / denorm2
        self.ny_mean_EMA = EMA(self.m * 1. / self.n)
        self.ny_stdev_EMA = EMA(self.m * 1. / self.n)
        
        self.start_tick = self.n * 4
        self.saliency = np.zeros(self.price.max() + self.simulator.stopprofit + 10)
        self.saliency2 = np.zeros(self.price.max() + self.simulator.stopprofit + 10)
        self.saliency_EMA = EMA(self.m)
        self.saliency2_EMA = EMA(self.m)
        self.mean_saliency_EMA = EMA(self.m * 1. / self.n)
        
        assert self.smooth_n % 2 == 1
        self.smooth_w = np.ones((self.smooth_n - 1) * self.tick_size + 1) / self.smooth_n 

def noiseFilter(data_in):

	N = int(np.ceil((4 / b)))
	if not N % 2: N += 1  # Make sure that N is odd.
	n = np.arange(N)
	 
	# Compute a low-pass filter with cutoff frequency fH.
	hlpf = np.sinc(2 * fH * (n - (N - 1) / 2.))
	hlpf *= np.blackman(N)
	hlpf = hlpf / np.sum(hlpf)
	 
	# Compute a high-pass filter with cutoff frequency fL.
	hhpf = np.sinc(2 * fL * (n - (N - 1) / 2.))
	hhpf *= np.blackman(N)
	hhpf = hhpf / np.sum(hhpf)
	hhpf = -hhpf
	hhpf[(N - 1) / 2] += 1
	 
	# Convolve both filters.
	h = np.convolve(hlpf, hhpf)
	s = np.convolve(data_in, hlpf)

	fig, ax = plt.subplots()
	ax.plot(data_in)
	plt.show()
	fig1, ax1 = plt.subplots()
	ax1.plot(s)
	plt.show()
	return s

def gauss_filter(dat, bin_freq, window=300, sigma=100):
    """
        turn psth into firing rate estimate. window size is in ms
    """
    if dat is None:
        return None, None
    window = np.int(1. / bin_freq * window)
    sigma = np.int(1. / bin_freq * sigma)
    r = range(-int(window / 2), int(window / 2) + 1)
    gaus = [1 / (sigma * np.sqrt(2 * np.pi)) *
            np.exp(-float(x) ** 2 / (2 * sigma ** 2)) for x in r]
    if len(dat.shape) > 1:
        fr = np.zeros_like(dat, dtype=np.float)
        for d in range(len(dat)):
            fr[d] = np.convolve(dat[d], gaus, 'same')
    else:
        fr = np.convolve(dat, gaus, 'same')
#    import pylab as plt
#    print bin_freq
#    plt.subplot(311)
#    plt.plot(gaus)
#    plt.subplot(312)
#    plt.plot(dat[:5].T)
#    plt.subplot(313)
#    plt.plot(fr[:5].T)
#    plt.show()

    return fr, len(gaus) / 2

def find_strand_shift(gdb, fwd_track, rev_track):
    # use a single chromosome to find shift between strands
    # that gives max covariance.
    if gdb.assembly == "dm3":
        chrom = gdb.get_chromosome("chr2L")
    else:
        chrom = gdb.get_chromosome("chr21")

    sys.stderr.write("retrieving values\n")
    fwd_vals = fwd_track.get_nparray(chrom)
    rev_vals = rev_track.get_nparray(chrom)

    # smooth values using sliding window
    sys.stderr.write("smoothing values\n")
    win = np.ones(SMOOTH_WIN_SIZE)
    fwd_vals = np.convolve(fwd_vals, win, mode="same")
    rev_vals = np.convolve(rev_vals, win, mode="same")

    # only use regions with high density of sites so all of contribution
    # to covariance comes from these
    fwd_cutoff = scipy.stats.scoreatpercentile(fwd_vals, 95)
    rev_cutoff = scipy.stats.scoreatpercentile(rev_vals, 95)
    fwd_vals[fwd_vals < fwd_cutoff] = 0
    rev_vals[rev_vals < rev_cutoff] = 0

    # find the shift that yields the max covariance
    max_cov_shift = find_max_cov(fwd_vals, rev_vals)

    return max_cov_shift

    def output(self):
        """The natural output for this analyzer is the analytic signal"""
        data = self.input.data
        sampling_rate = self.input.sampling_rate
        
        a_signal =\
    ts.TimeSeries(data=np.zeros(self.freqs.shape+data.shape,
                                        dtype='D'),sampling_rate=sampling_rate)
        if self.freqs.ndim == 0:
            w = self.wavelet(self.freqs,self.sd,
                             sampling_rate=sampling_rate,ns=5,
                                                     normed='area')

            nd = (w.shape[0]-1)/2
            a_signal.data[...] = (np.convolve(data,np.real(w),mode='same') +
                                  1j*np.convolve(data,np.imag(w),mode='same'))
        else:    
            for i,(f,sd) in enumerate(zip(self.freqs,self.sd)):
                w = self.wavelet(f,sd,sampling_rate=sampling_rate,
                                 ns=5,normed='area')

                nd = (w.shape[0]-1)/2
                a_signal.data[i,...] = (np.convolve(data,np.real(w),mode='same')+1j*np.convolve(data,np.imag(w),mode='same'))
                
        return a_signal

 def apply_window_function(cls, input_array, window_size, output_array):
     sum_filter = np.ones((window_size,), dtype=np.float32)/window_size
     squares = np.square(input_array)
     sum_of_squares = np.convolve(squares, sum_filter, mode='valid')
     sums = np.convolve(input_array, sum_filter, mode='valid')
     squares_of_sums = np.square(sums)
     output_array[:] = squares_of_sums/sum_of_squares

def slidingWindowV(P,inner=3,outer=64,maxM=50,minM=7,maxT=59,norm=True):
	""" Enhance the constrast vertically (along frequency dimension)

		Cut off extreme values and demean the image
		Utilize numpy convolve to get the mean at a given pixel
		Remove local mean with inner exclusion region

		Args:
			P: 2-d numpy array image
			inner: inner exclusion region 
			outer: length of the window
			maxM: size of the output image in the y-dimension
			norm: boolean to cut off extreme values

		Returns:
			Q: 2-d numpy contrast enhanced vertically
	"""
	Q = P.copy()
	m, n = Q.shape
		
	if norm:
		mval, sval = np.mean(Q[minM:maxM,:maxT]), np.std(Q[minM:maxM,:maxT])
		fact_ = 1.5
		Q[Q > mval + fact_*sval] = mval + fact_*sval
		Q[Q < mval - fact_*sval] = mval - fact_*sval
		Q[:minM,:] = mval
	wInner = np.ones(inner)
	wOuter = np.ones(outer)
	for i in range(maxT):
		Q[:,i] = Q[:,i] - (np.convolve(Q[:,i],wOuter,'same') - np.convolve(Q[:,i],wInner,'same'))/(outer - inner)
	Q[Q < 0] = 0.

	return Q[:maxM,:]

def get_spline(data):
	""" Returns array of cubic spline interpolation polynoms (for every interval [x[i -1]; x[i]] ) """
	h = set_h(data)
	A = set_matrix_A(h)
	B = set_vector_B(data, h)
	m = find_vector_m(A, B)	
	spline_array = []

	for i in range(1,len(data)):
		xi_minus_x_cub = [(data[i][0] ** 3), -3 * (data[i][0] ** 2), 3*(data[i][0]),   -1]
		s1 = list(np.convolve( (m[i - 1]  / (6 * h[i-1])), xi_minus_x_cub))	
		x_minus_xi_1_cub = [-(data[i-1][0] ** 3), 3 * (data[i-1][0] ** 2), -3 * data[i-1][0],  1]
		s2 = list(np.convolve( (m[i] / (6 * h[i-1])), x_minus_xi_1_cub ))
		ai = data[i-1][1] - ((m[i-1]*h[i-1]**2)/6)
		s3 = list(np.convolve((ai/h[i-1]), [data[i][0], -1]))
		bi = data[i][1] - ((m[i]*h[i-1]**2)/6)
		s4 = list(np.convolve((bi/h[i-1]), [-data[i-1][0], 1]))

		iter_length = max(len(s1), len(s2), len(s3), len(s4))
		for k in range(iter_length - len(s1)): s1.append(0)
		for k in range(iter_length - len(s2)): s2.append(0)
		for k in range(iter_length - len(s3)): s3.append(0)
		for k in range(iter_length - len(s4)): s4.append(0)

		spline = [0 for t in range(iter_length)]
		for j in range(iter_length):
			spline[j] = s1[j] + s2[j] + s3[j] + s4[j] 
		spline_array.append(spline)
	return spline_array

def perbaseToWindowAverageMasked( 
                            V, 
                            arMasked,
                            cooStart,
                            cooStop,
                            wndWidth,
                            kernel='sum',  # can be sum, or mean
                           ):
    assert arMasked.shape==V.shape, (arMasked.shape,V.shape)
    
    
    krn = np.ones( (wndWidth,), 'int32' )
        
    Vconv = np.convolve( V*((1-arMasked).astype('uint64')), krn, 'full' )
    Vconv = Vconv[ (wndWidth-1):Vconv.shape[0]:wndWidth ]
        
    assert Vconv.shape[0]==(V.shape[0]/wndWidth + min(V.shape[0]%wndWidth,1))
        
    nBpMasked = np.convolve( arMasked, np.ones(wndWidth,'int32'), 'full' )
    nBpMasked = nBpMasked[ (wndWidth-1):nBpMasked.shape[0]:wndWidth ]
    
    assert nBpMasked.shape[0]==Vconv.shape[0]
    
    outStarts = cooStart+np.arange(Vconv.shape[0])*wndWidth
    
    nBpUnmasked = wndWidth - nBpMasked

    if kernel=='mean':
        nBpUnmasked[-1] = cooStop - outStarts[-1] + 1
        
        Vconv = (ma.masked_array( Vconv, nBpUnmasked == 0 ) / nBpUnmasked.astype('float64')).filled(0.)
        
    return outStarts,Vconv,nBpUnmasked

def score_region(motifs, seq):
    # code the sequence
    coded_seq = {}
    for base in 'ACGT':
        coded_seq[base] = np.zeros(len(seq), dtype='float32')
    for i, base in enumerate(seq.upper()):
        coded_seq[base][i] = 1
    coded_RC_seq = {}
    for base in 'TGCA':
        coded_seq[base] = np.zeros(len(seq), dtype='float32')


    motif_scores = []
    for motif in motifs:
        score_mat = motif.motif_data
        scores = np.zeros(len(seq)-len(score_mat)+1, dtype='float32')
        for base, base_scores in zip('ACGT', score_mat.T):
            scores += np.convolve(coded_seq[base], base_scores, mode='valid')

        for i, base in enumerate(seq.upper()):
            coded_seq[base][len(seq)-i-1] = 1
        RC_scores = np.zeros(len(seq)-len(score_mat)+1, dtype='float32')
        for base, base_scores in zip('TCGA', score_mat.T):
            scores += np.convolve(coded_seq[base], base_scores, mode='valid')

        max_scores = np.vstack((scores, RC_scores)).max(0)
        motif_scores.append( max_scores.mean() )

    return np.array(motif_scores)

    def richardson_lucy_deconvolution(self,  psf, iterations=15, 
                                      mask=None):
        """1D Richardson-Lucy Poissonian deconvolution of 
        the spectrum by the given kernel.
    
        Parameters
        ----------
        iterations: int
            Number of iterations of the deconvolution. Note that 
            increasing the value will increase the noise amplification.
        psf: EELSSpectrum
            It must have the same signal dimension as the current 
            spectrum and a spatial dimension of 0 or the same as the 
            current spectrum.
            
        Notes:
        -----
        For details on the algorithm see Gloter, A., A. Douiri, 
        M. Tence, and C. Colliex. “Improving Energy Resolution of 
        EELS Spectra: An Alternative to the Monochromator Solution.” 
        Ultramicroscopy 96, no. 3–4 (September 2003): 385–400.
        
        """
        self._check_signal_dimension_equals_one()
        ds = self.deepcopy()
        ds.data = ds.data.copy()
        ds.mapped_parameters.title += (
            ' after Richardson-Lucy deconvolution %i iterations' % 
                iterations)
        if ds.tmp_parameters.has_item('filename'):
                ds.tmp_parameters.filename += (
                    '_after_R-L_deconvolution_%iiter' % iterations)
        psf_size = psf.axes_manager.signal_axes[0].size
        kernel = psf()
        imax = kernel.argmax()
        j = 0
        maxval = self.axes_manager.navigation_size
        if maxval > 0:
            pbar = progressbar(maxval=maxval)
        for D in self:
            D = D.data.copy()
            if psf.axes_manager.navigation_dimension != 0:
                kernel = psf(axes_manager=self.axes_manager)
                imax = kernel.argmax()

            s = ds(axes_manager=self.axes_manager)
            mimax = psf_size -1 - imax
            O = D.copy()
            for i in xrange(iterations):
                first = np.convolve(kernel, O)[imax: imax + psf_size]
                O = O * (np.convolve(kernel[::-1], 
                         D / first)[mimax: mimax + psf_size])
            s[:] = O
            j += 1
            if maxval > 0:
                pbar.update(j)
        if maxval > 0:
            pbar.finish()
        
        return ds

def high_pass_filter(x, y, L):
    y_filter = -np.ones(L)/L
    y_filter[L/2] += 1
    new_y = np.convolve(y, y_filter, 'valid') 
    x_filter = np.zeros(L)
    x_filter[L/2] = 1
    return np.convolve(x, x_filter, 'valid'), np.convolve(y, y_filter, 'valid')

def running_average_masked(obs, ws, min_valid_fraction=0.95):
    '''
    calculates a running average via convolution, fixing the edges
    obs     --  observations (a masked array)
    ws      --  window size (number of points to average)
    '''
    #tmp_vals = np.convolve(np.ones(ws, dtype=float), obs*(1-obs.mask), mode='same')
    tmp_vals = np.convolve(np.ones(ws, dtype=float), obs.filled(0), mode='same')

     # if the array is not masked, edges needs to be explictly fixed due to smaller counts
    if len(obs.mask.shape) == 0:
        tmp_valid = ws*np.ones_like(tmp_vals)
        # fix the edges. using mode='same' assumes zeros outside the range
        if ws%2==0:
            tmp_vals[:ws//2]*=float(ws)/np.arange(ws//2,ws)
            if ws//2>1:
                tmp_vals[-ws//2+1:]*=float(ws)/np.arange(ws-1,ws//2,-1.0)
        else:
            tmp_vals[:ws//2]*=float(ws)/np.arange(ws//2+1,ws)
            tmp_vals[-ws//2:]*=float(ws)/np.arange(ws,ws//2,-1.0)

    # if the array is masked, then we get the normalizer from counting the unmasked values
    else:
        tmp_valid = np.convolve(np.ones(ws, dtype=float), (1-obs.mask), mode='same')

    run_avg = np.ma.array(tmp_vals / tmp_valid)
    run_avg.mask = tmp_valid < ws * min_valid_fraction

    return run_avg

def plot_data(data, nplots = 4):
    global window_size_
    window = np.ones(window_size_) / window_size_
    tvec, yvec = data[:,0], data[:,1]
    pylab.subplot(nplots, 1, 1)
    pylab.plot(tvec, yvec, label="raw data")
    pylab.legend()

    yvec = np.convolve(yvec, window, 'same')
    pylab.subplot(nplots, 1, 2)
    pylab.plot(tvec, yvec, label='Window size = %s' % window_size_)
    pylab.plot([0, tvec[-1]], [0.5*np.mean(yvec)]*2, label = '0.5*Mean pupil size')
    pylab.legend()

    pylab.subplot(nplots, 1, 4)
    # When area reduces to half of eye pupil, it should be considered.
    newY = 0.5*yvec.mean() - yvec
    newY = newY + np.fabs(newY)
    window = np.ones(3*window_size_)/(3*window_size_)

    yy = np.convolve(newY, window, 'same')
    pylab.plot(tvec, yy, label='Blinks')

    pylab.xlabel("Time (seconds)")
    outfile = 'output.png'
    print("[INFO] Writing to %s" % outfile)
    pylab.savefig(outfile)

def csd():
    import numpy as np
    from smartplotlib import xyplot, alias
    # make some data
    dt = 0.01
    t = np.arange(0, 30, dt)
    nse1 = np.random.randn(len(t))                 # white noise 1
    nse2 = np.random.randn(len(t))                 # white noise 2
    r = np.exp(-t/0.05)

    cnse1 = np.convolve(nse1, r, mode='same')*dt   # colored noise 1
    cnse2 = np.convolve(nse2, r, mode='same')*dt   # colored noise 2

    # two signals with a coherent part and a random part
    s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1
    s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2

    xy = xyplot.derive(xlabel="time", ylabel="s1 & s2",
                       xlim=(0,5), axes=211)


    csd = xyplot.csd(s1, s2, 256, 1./dt, axes=212)

    xy.go("fclear", "axes")
    xy.plot(t, s1, 'b-', t, s2, 'g-')
    csd.go("plot", "axes", "show", "draw")

def find_blinks_using_edge(data, plot = False, **kwargs):
    """Find location of blinks in data"""
    global window_size_
    records = OrderedDict()
    window = np.ones(window_size_)/window_size_
    t, y = data[:,0], data[:,1]
    # Smooth out the vectors.
    yvec = np.convolve(y, window, 'same')
    records['smooth'] = (t, y)
    newY = 0.5*yvec.mean() - yvec
    newY = newY + np.fabs(newY)
    window = np.ones(window_size_)/(window_size_)
    yy = np.convolve(newY, window, 'same')
    blinks = []
    while yy.max() > 10:
        i = np.argmax(yy)
        isBlink, a = get_blink(i, yy)
        if isBlink:
            blinks.append((i, a))

    xvec, yvec = [], []
    for i, x in sorted(blinks):
        xvec.append(t[i])
        yvec.append(x)
    return xvec, yvec

def cohere():
    """
    Compute the coherence of two signals
    mpl_examples/pylab_examples/cohere_demo.py
    """
    import numpy as np
    from smartplotlib import xyplot, subplot, alias


    # make a little extra space between the subplots

    dt = 0.01
    t = np.arange(0, 30, dt)
    nse1 = np.random.randn(len(t))                 # white noise 1
    nse2 = np.random.randn(len(t))                 # white noise 2
    r = np.exp(-t/0.05)

    cnse1 = np.convolve(nse1, r, mode='same')*dt   # colored noise 1
    cnse2 = np.convolve(nse2, r, mode='same')*dt   # colored noise 2

    # two signals with a coherent part and a random part
    s1 = 0.01*np.sin(2*np.pi*10*t) + cnse1
    s2 = 0.01*np.sin(2*np.pi*10*t) + cnse2

    ps1 = xyplot(t, s1, fmt='b-', y2= s2, fmt2="g-",
                 axes=211, xlabel="time", ylabel="s1 & s2")

    ps1.go("fclear", "axes", "plot")

    c = xyplot.cohere(s1, s2, 256, 1./dt,
                      ylabel="coherence", axes=212)

    c.go("axes", "plot", "draw", "show")
    return c

 def smooth(self, signal, pad=True):
     """
     Returns smoothed signal (or it's n-th derivative).
     
     Parameters
     ---------- 
     y : array_like, shape (N,)
         the values of the time history of the signal.
     pad : bool
        pad first and last values to lessen the end effects.
         
     Returns
     -------
     ys : ndarray, shape (N)
         the smoothed signal (or it's n-th derivative).
     """
     coeff = self._coeff
     n = size(coeff - 1) // 2
     y = np.squeeze(signal)
     if pad:
         first_vals = y[0] - abs(y[n:0:-1] - y[0])
         last_vals = y[-1] + abs(y[-2:-n - 2:-1] - y[-1]) 
         y = concatenate((first_vals, y, last_vals))
         n *= 2 
     d = y.ndim
     if d > 1:
         y1 = y.reshape(y.shape[0], -1)
         res = []
         for i in range(y1.shape[1]):
             res.append(convolve(y1[:, i], coeff)[n:-n])
         res = np.asarray(res).T
     else:
         res = convolve(y, coeff)[n:-n]
     return res

 def _jacobian(self, param, y, weights=None):
     if weights is None:
         weights = 1.
     if self.convolved is True:
         counter = 0
         grad = np.zeros(len(self.axis.axis))
         for component in self:  # Cut the parameters list
             if component.active:
                 component.fetch_values_from_array(
                     param[
                         counter:counter +
                         component._nfree_param],
                     onlyfree=True)
                 if component.convolved:
                     for parameter in component.free_parameters:
                         par_grad = np.convolve(
                             parameter.grad(self.convolution_axis),
                             self.low_loss(self.axes_manager),
                             mode="valid")
                         if parameter._twins:
                             for par in parameter._twins:
                                 np.add(par_grad, np.convolve(
                                     par.grad(
                                         self.convolution_axis),
                                     self.low_loss(self.axes_manager),
                                     mode="valid"), par_grad)
                         grad = np.vstack((grad, par_grad))
                 else:
                     for parameter in component.free_parameters:
                         par_grad = parameter.grad(self.axis.axis)
                         if parameter._twins:
                             for par in parameter._twins:
                                 np.add(par_grad, par.grad(
                                     self.axis.axis), par_grad)
                         grad = np.vstack((grad, par_grad))
                 counter += component._nfree_param
         to_return = grad[1:, self.channel_switches] * weights
     else:
         axis = self.axis.axis[self.channel_switches]
         counter = 0
         grad = axis
         for component in self:  # Cut the parameters list
             if component.active:
                 component.fetch_values_from_array(
                     param[
                         counter:counter +
                         component._nfree_param],
                     onlyfree=True)
                 for parameter in component.free_parameters:
                     par_grad = parameter.grad(axis)
                     if parameter._twins:
                         for par in parameter._twins:
                             np.add(par_grad, par.grad(
                                 axis), par_grad)
                     grad = np.vstack((grad, par_grad))
                 counter += component._nfree_param
         to_return = grad[1:, :] * weights
     if self.signal.metadata.Signal.binned is True:
         to_return *= self.signal.axes_manager[-1].scale
     return to_return