Python numpy.convolve函数代码示例

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)