def vsumsmooth_offset(x, w, center, csum, delta, offset,spacing = 0.01):
	exact=vsumexact_offset(x,w,center,csum,offset)
	if np.shape(x)==():
		gaussian = signal.gaussian(10*delta/spacing,delta/spacing)/(delta/spacing*np.sqrt(2*np.pi))
	else:
		gaussian = signal.gaussian(10*delta/(x.values[1]-x[0]),delta/((x.values[1]-x[0])))/(delta/(x.values[1]-x[0])*np.sqrt(2*np.pi))
	return signal.fftconvolve(exact,gaussian,mode='same')

def vdifsmooth(x, w, center, cdif, delta,spacing=0.01):
	exact=vdifexact(x,w,center,cdif)
	if np.shape(x) == ():
		gaussian = signal.gaussian(10*delta/spacing,delta/spacing)/(delta/spacing*np.sqrt(2*np.pi))
	else:
		gaussian = signal.gaussian(10*delta/(x.values[1]-x[0]),delta/((x.values[1]-x[0])))/(delta/(x.values[1]-x[0])*np.sqrt(2*np.pi))
	return signal.fftconvolve(exact,gaussian,mode='same')

def makescore(temp,dims):
    score = zeros(dims) + 1
    
    lonax = greater(dims[1],dims[0])
    londim= dims[lonax]
    sdim = dims[1-lonax]
    ldiag = arange(0,londim,1,int)
    sdiag = array(floor(ldiag * (float(sdim)/londim)),int)
    dvals = [ldiag, sdiag]
    if lonax:dvals = dvals[::-1]
    score[dvals] = 20.*temp ** .5

    if lonax: score = score.T

    #for i in range(len(score)):
    #    score[i][i] = 20.* temp **.5

    if temp > .1:
        g = ss.gaussian((londim/2)*temp**2,(londim/2)*temp**2)[:,newaxis]*\
            ss.gaussian((sdim/2)*temp**2,(sdim/2)*temp**2)[newaxis,:]

        g/= sum(g)
        score = ss.convolve2d(score,g,'same','wrap')
        
    for i in range(1,len(score)):
        score[i,arange(sdiag[i])] *= .25


        
    if lonax: score = score.T

    return score

def fftgauss(img,sigma):

    """https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.signal.fftconvolve.html"""

    kernel = np.outer(signal.gaussian(img.shape[0], sigma),
                        signal.gaussian(img.shape[1],sigma))

    return signal.fftconvolve(img, kernel, mode='same')

def addGaussians(query_fn,cand_list,tree,K):
	songs_list   = getNeighbors(query_fn,cand_list,tree,K)
	M            = 23000
	query_ann    = getAnnotationList(gt_path,[query_fn])
	query_labels = [elem[-1] for elem in query_ann[0]]
	query_ann    = np.floor((np.array(getAnnotation(query_ann))*1000)).astype(int)
	length       = query_ann[-1]
	total        = np.zeros(int(np.ceil(length)))

	neighbors_annotations_rescaled = []
	neighbors_annotations = getAnnotationList(gt_path,songs_list)

	for i, song in enumerate(songs_list):
		gt_list        = getAnnotationList(gt_path,[song])
		ann            = np.floor((np.array(getAnnotation(gt_list))*1000)).astype(int) #convert to miliseconds to mantain res
		neighbor_dur   = ann[-1]
		ann_with_sides = ann
		ann            = ann[1:-1]
		a              = np.zeros(int(np.ceil(length)))
		r              = float(length)/float(neighbor_dur) #rescale according to query duration
		ann            = np.floor(ann*r)
		ann_with_sides = np.floor(ann_with_sides*r) 

		labels = [x[-1] for x in gt_list[0]] # get the labels
		annotation_rescaled = []
		for elem in neighbors_annotations[i]:
			label = elem[-1] #save the label so it doesnt get affected by rescaling
			elem[0] = int(np.floor(float(elem[0])*1000*r)) #rescale the rest
			elem[1] = int(np.floor(float(elem[1])*1000*r))
			annotation_rescaled.append([elem[0],elem[1],label])
		neighbors_annotations_rescaled.append(annotation_rescaled)
		for i, loc in enumerate(ann,1):
			section_length = ann_with_sides[i]-ann_with_sides[i-1]
			sigma = 0.1*section_length
			# M=int(np.floor(0.6*section_length))
			g1 = signal.gaussian(M,std=sigma)
			half1 = int(np.floor(len(g1)/2))
			section_length = ann_with_sides[i+1]-ann_with_sides[i]
			sigma = 0.1*section_length
			g2 = signal.gaussian(M,std=sigma)
			half2 = int(np.floor(len(g2)/2))
			g = np.concatenate((g1[:half1],g2[half2:]))			
			if loc < np.floor(M/2):
				a += np.array(np.concatenate((g[int(np.floor(M/2)-loc):],np.zeros(int(length-loc-np.floor(M/2))))))
			elif loc + np.floor(M/2) > length:
				a += np.array(np.concatenate((np.zeros(int(loc-np.floor(M/2))),g[:int(length+np.floor(M/2)-loc)])))
			else:
				a += np.array(np.concatenate((np.zeros(int(loc-np.floor(M/2))),g,np.zeros(int(length-loc-np.floor(M/2))))))
		total += a
	total = total/float(max(total))
	peaks = getPeaks(total,neighbors_annotations)
	all_songs_segmented = [segmentLabel(elem) for elem in neighbors_annotations_rescaled]
	res_boundaries = sorted(peaks)
	res_boundaries.insert(0,0)
	res_boundaries.append(length)
	res_labels = mergeLabels(res_boundaries,all_songs_segmented)
	res_annotations = formatAnnotation(res_boundaries,res_labels)
	return res_annotations

    def _gaussian_window(self, width, sigma):
        """
        Generates a gaussian window

        sigma is based on the dat being in a range 0 to 1
        """
        return (ssignal.gaussian(width[0], sigma*width[0]).reshape((-1,1,1)) *
                ssignal.gaussian(width[1], sigma*width[1]).reshape((-1,1)) *
                ssignal.gaussian(width[2], sigma*width[2]))

def heat_labels_gauss(click, img_size=IMG_SIZE, k_size=KERNEL_SIZE, label_size=LABEL_SIZE):
    # take list of pixel coordinates and return 70x70 heatmap
    img = np.zeros((img_size, img_size))
    for j in range(click.shape[0]):
        x = img_size-1-click[j,1]
        y = click[j,0]
        img[x,y]=1
    kernel = np.outer(signal.gaussian(img_size+1, k_size), signal.gaussian(img_size+1, k_size))
    img = signal.convolve2d(img,kernel, mode='same')
    offset = (img_size-img_size/label_size*(label_size-1))/2
    step = img_size/label_size
    return img[offset:(img_size-offset+step):step, offset:(img_size-offset+step):step]

def smooth_pdf(a, sd=None):
    """Get a smoothed pdf of an array of data for visualization
    
    Keyword arguments:
    sd -- S.D. of the gaussian kernel used to perform the smoothing (default is
        1/20 of the data range)
    
    Return 2-row (x, pdf(x)) smoothed probability density estimate.
    """
    from scipy.signal import gaussian, convolve
    from numpy import array, arange, cumsum, trapz, histogram, diff, r_, c_
    if sd is None:
        sd = 0.05 * a.ptp()
    data = a.copy().flatten() # get 1D copy of array data
    nbins = len(data) > 1000 and len(data) or 1000 # num bins >~ O(len(data))
    f, l = histogram(data, bins=nbins, normed=True) # fine pdf
    sd_bins = sd * (float(nbins) / a.ptp()) # convert sd to bin units
    kern_size = int(10*sd_bins) # sufficient convolution kernel size
    g = gaussian(kern_size, sd_bins) # generate smoothing kernel
    c = cumsum(f, dtype='d') # raw fine-grained cdf
    cext = r_[array((0,)*(2*kern_size), 'd'), c, 
        array((c[-1],)*(2*kern_size), 'd')] # wrap data to kill boundary effect
    cs = convolve(cext, g, mode='same') # smooth the extended cdf
    ps = diff(cs) # differentiate smooth cdf to get smooth pdf
    dl = l[1]-l[0] # get bin delta
    l = r_[arange(l[0]-kern_size*dl, l[0], dl), l, 
        arange(l[-1]+dl, l[-1]+kern_size*dl, dl)] # pad index to match bounds
    ps = ps[kern_size:kern_size+len(l)] # crop pdf to same length as index
    ps /= trapz(ps, x=l) # normalize pdf integral to unity
    return c_[l, ps].T # return 2-row concatenation of x and pdf(x)

def smooth_color_prior(size=64, sigma=5, do_plot=False):

    prior_prob = np.load(os.path.join(data_dir, "CelebA_%s_prior_prob.npy" % size))
    # add an epsilon to prior prob to avoid 0 vakues and possible NaN
    prior_prob += 1E-3 * np.min(prior_prob)
    # renormalize
    prior_prob = prior_prob / (1.0 * np.sum(prior_prob))

    # Smooth with gaussian
    f = interp1d(np.arange(prior_prob.shape[0]),prior_prob)
    xx = np.linspace(0,prior_prob.shape[0] - 1, 1000)
    yy = f(xx)
    window = gaussian(2000, sigma)  # 2000 pts in the window, sigma=5
    smoothed = convolve(yy, window / window.sum(), mode='same')
    fout = interp1d(xx,smoothed)
    prior_prob_smoothed = np.array([fout(i) for i in range(prior_prob.shape[0])])
    prior_prob_smoothed = prior_prob_smoothed / np.sum(prior_prob_smoothed)

    # Save
    file_name = os.path.join(data_dir, "CelebA_%s_prior_prob_smoothed.npy" % size)
    np.save(file_name, prior_prob_smoothed)

    if do_plot:
        plt.plot(prior_prob)
        plt.plot(prior_prob_smoothed, "g--")
        plt.plot(xx, smoothed, "r-")
        plt.yscale("log")
        plt.show()

def make_gaussian(k, std):
  '''Create a gaussian kernel.
  
  Input:

  k - the radius of the kernel.
  
  std - the standard deviation of the kernel.
  
  Output:

  output - a numpy array of shape (2k+1, 2k+1) and dtype float.
  
  If gaussian_1d is a gaussian filter of length 2k+1 in one dimension, 
  kernel[i,j] should be filled with the product of gaussian_1d[i] and 
  gaussian_1d[j].
 
  Once all the points are filled, the kernel should be scaled so that the sum
  of all cells is equal to one.'''
  kernel = np.zeros((2*k+1, 2*k+1), dtype = float)
  gaussian1d = signal.gaussian(2*k+1, std)
  for i in range(0,kernel.shape[0]):
    for j in range(0,kernel.shape[1]):
       kernel[i,j] = gaussian1d[i]*gaussian1d[j]
  kernel = kernel/kernel.sum()

  # Insert your code here.----------------------------------------------------

  #---------------------------------------------------------------------------
  return kernel

def process_dir(sDir, iResample=None, iSmooth = 50, iSigmaSecs=0.01):
    """
    take input dir and output smoothed, correlated array
    iSigmaSecs:  standard deviation of gaussian in seconds
    iSmooth = smoothing window size for linear smoother
    """
    from scipy import signal
    import numpy as np
    iSampleRate, aTime, aOrigAudio = audio2array(sDir, iResample)
    
    #only positive
    aAudio = [abs(i) for i in aOrigAudio]
    
    #audio files must be right format
    aOrigAudio = np.asarray(aOrigAudio, dtype=np.int16)
    
    if not iSmooth == None:
        #smooth
        aAudio = smooth(aAudio, iSmooth)
    
    #standard deviation for gaussian function
    iSigma = float(iSigmaSecs * iSampleRate)
    aGaussian = signal.gaussian(10*iSigma, iSigma)
    
    #gaussian correlated with audio signal
    aCorr = np.correlate(aAudio, aGaussian, 'same')
    
    
    return iSampleRate, aTime, aAudio, aCorr, aOrigAudio

 def simulSpectra(self,fileSpectra, nBin, continuumLevel, sigma, linePositionChan, lineWidth, lineIntensity, startFreq = 0., resFreq =1.):
     """ Save in fileSpectra of nBin channels with the linePosition, lineWidth and lineIntensity list 
     starFreq and resFreq are optional"""
     
     freq = np.arange(nBin)
     spectra = np.random.normal(continuumLevel, sigma, nBin)
     
     index = 0
     for pos in linePositionChan:
         
         nChan = 4 * int(lineWidth[index] / resFreq)    
         spec = lineIntensity[index] * signal.gaussian(nChan,lineWidth[index])
         
         startPos = pos - nChan / 4
         spectra[pos:pos+nChan] = spectra[pos:pos+nChan] + spec
     
         index += 1
         
         
     f = open(fileSpectra,"w")
     
     index = 0
     for frequency  in freq :
         strOut = "%f   %f \n"%(frequency, spectra[index])
         f.write(strOut)
         index += 1
         
     f.close()

def make_gaussian(k, std):
  '''Create a gaussian kernel.
  
  Input:

  k - the radius of the kernel.
  
  std - the standard deviation of the kernel.
  
  Output:

  output - a numpy array of shape (2k+1, 2k+1) and dtype float.
  
  If gaussian_1d is a gaussian filter of length 2k+1 in one dimension, 
  kernel[i,j] should be filled with the product of gaussian_1d[i] and 
  gaussian_1d[j].
 
  Once all the points are filled, the kernel should be scaled so that the sum
  of all cells is equal to one.'''
  kernel = None
  # Insert your code here.----------------------------------------------------
  l = 2 * k + 1
  kernel = np.zeros(l * l)
  kernel.shape = (l, l)
  from scipy import signal
  gaussian1d = signal.gaussian(l, std)
  for i in xrange(l):
    for j in xrange(l):
      kernel[i, j] = gaussian1d[i] * gaussian1d[j]
  s = np.sum(kernel)
  kernel = np.divide(kernel, s)
  #---------------------------------------------------------------------------
  return kernel

def moving_average(series,sigma = 3,window_time = 39):
    #### Moving weighted gaussian average with window = 39
    b = gaussian(window_time,sigma)
    average = filters.convolve1d(series,b/b.sum())
    var = filters.convolve1d(np.power(series-average,2),b/b.sum())
    
    return average,var

def compute_gaussian_krnl(M):
    """Creates a gaussian kernel following Foote's paper."""
    g = signal.gaussian(M, M / 3., sym=True)
    G = np.dot(g.reshape(-1, 1), g.reshape(1, -1))
    G[M / 2:, :M / 2] = -G[M / 2:, :M / 2]
    G[:M / 2, M / 2:] = -G[:M / 2, M / 2:]
    return G

def testGauss(x, y, s, npts):
	#b = gaussian(39, 10)
	b = gaussian(75, 15)
	ga = filters.convolve1d(y, b/b.sum())
	plt.plot(x, ga)
	print "gaerr", ssqe(ga, s, npts)
	return ga

    def __init__(self, brain):
        """A convergence-divergence zone.
        1. correlates temporally near packets and store the correlations allowing for quick and efficient lookups.
        2. Allows for generating of 2nd modality output given first modality.
        """
        self.brain = brain
        self.LEARNING_RATE = config.CE_LEARNING_RATE

        # The maximum number of timesteps we look back.
        # IE. The point where we trim the tails of the Gaussian
        # This is for computational efficiency
        gaussian = signal.gaussian(config.CE_CORRELATION_WINDOW_MAX * 2, std=config.CE_CORRELATION_WINDOW_STD, sym=True)
        self.GAUSSIAN = np.split(gaussian, 2)[0][::-1]  # split the array into two and then reverse it
        assert len(self.GAUSSIAN) == config.CE_CORRELATION_WINDOW_MAX

        if config.CE_IGNORE_GAUSSIAN:
            self.GAUSSIAN *= 0
            self.GAUSSIAN[0] = 1

        # The number of packets to keep in the queue
        # The number of packets used for learning from the packet_queue depends on the
        self.PACKET_QUEUE_LENGTH = len(self.GAUSSIAN) + 1
        # ================== END CONFIG ==================================

        # A queue to store the most recent packets, old packets are automatically pushed out of the queue.
        self.packet_queue = deque(maxlen=self.PACKET_QUEUE_LENGTH)

        # A dict to store the connections/correlations between the different modalities
        self.correlations = {}

def gbm(R,sigma):
    # gaussian_blur_matrix
    #   R   为模糊半径
    # sigma 为标准差
    temp1=signal.gaussian(2*R+1,sigma)
    temp2=np.sum(temp1)  #用于归一化，使2R+1矩阵的和为1
    return temp1/temp2

def spectrogram_scores(pattern_chunk, chan_sound, candidates):
    s_f = chan_sound.s_f
    n_window = 256
    n_overlap = 192
    sigma = 1. / 1000. * s_f

    spectrogram_kwargs = {'nperseg': n_window,
                          'noverlap': n_overlap,
                          'window': sg.gaussian(n_window, sigma),
                          'scaling': 'spectrum'}

    pattern_spectrogram = spectrogram(pattern_chunk.data[:, 0], s_f, **spectrogram_kwargs)

    logger.info('Getting spectrogram difference score for {} candidates'.format(len(candidates.index)))
    for (i, start) in enumerate(candidates['start'][:]):
        logger.debug('Start {0}: {1}'.format(i, start))
        motif_start = start
        series = chan_sound.get_chunk(motif_start, motif_start + pattern_chunk.samples)
        # f, t, sxx = spectrogram(bandpass_filter(series[:, 0], s_f), s_f, **spectrogram_kwargs)

        candidates.set_value(i, 'spectral_diff',
                             spectrogram_diff(series[:, 0],
                                              pattern_spectrogram[2],
                                              s_f,
                                              spectrogram_kwargs)
                             )

def make_gaussian(k, std):
  '''Create a gaussian kernel.

  Input:

  k - the radius of the kernel.

  std - the standard deviation of the kernel.

  Output:

  output - a numpy array of shape (2k+1, 2k+1) and dtype float.

  If gaussian_1d is a gaussian filter of length 2k+1 in one dimension,
  kernel[i,j] should be filled with the product of gaussian_1d[i] and
  gaussian_1d[j].

  Once all the points are filled, the kernel should be scaled so that the sum
  of all cells is equal to one.'''
  kernel = None
  # Insert your code here.----------------------------------------------------
  size = 2 * k + 1
  gaussian1d = signal.gaussian(size, std)
  gaussian2d = np.ndarray((size, 1), buffer=gaussian1d) * np.ndarray((1, size), buffer=gaussian1d)
  kernel = gaussian2d / np.sum(gaussian2d)
  #---------------------------------------------------------------------------
  return kernel