def degrades(V, W, H, rate, s, niter=20):
    """
    Deconvolution by Gradient Descent with Sparsification.
    """
    V, W, H = normalize(V), normalize(W), normalize(H)
    for iter in xrange(niter):
        convolved_pieces = np.vstack([signal.convolve(W[r], H[r], mode='full') for r in xrange(W.shape[0])])
        convolved = np.sum(convolved_pieces, axis=0)
        delta = V - convolved
        projected_H = fft_deconvolve(delta, W)
        H = magical_entropy_hammer(H, H + projected_H * rate, 0.01, axis=1)
        
        
        convolved_pieces = np.vstack([signal.convolve(W[r], H[r], mode='full') for r in xrange(W.shape[0])])
        convolved = np.sum(convolved_pieces, axis=0)
        delta = V - convolved
        projected_W = fft_deconvolve(delta, H)
        W = magical_entropy_hammer(W, W + projected_W * rate, 0.001, axis=1)
        
        print np.sum(np.abs(delta))
        print H
        
    pylab.clf()
    pylab.subplot(311)
    pylab.plot(W.T)
    pylab.subplot(312)
    pylab.plot(H.T)
    pylab.subplot(313)
    pylab.plot(V)
    pylab.plot(convolved)
    return W, H

def test3DConvolution():
    kernel = array([[[-1, 1], [1, 0]], [[1, 0], [0, 0]]])
    kernel_cube = array([
        zeros((3,3)),
        [[0, 0, 0], [0, -1, 1], [0, 1, 0]],
        [[0, 0, 0], [0,  1, 0], [0, 0, 0]]])

    A = reshape(linspace(1,9,9), (3,3), order = 'C')
    B = array((A, A, A))

    C = operators.convolve(kernel, (3, 3, 3), order = 'F')
    adjointTest(C)

    np.testing.assert_allclose(
        reshape(C._forward(ravel(B, order = 'F')), (3, 3, 3), order = 'F'),
        signal.convolve(B, kernel, 'same'))

    np.testing.assert_allclose(
        reshape(C._forward(ravel(B, order = 'F')), (3, 3, 3), order = 'F'),
        signal.convolve(B, kernel_cube, 'same'))

    np.testing.assert_allclose(
        reshape(C._adjoint(ravel(B, order = 'F')), (3, 3, 3), order = 'F'),
        signal.convolve(B,
            np.flipud(np.fliplr(kernel_cube[:,:,::-1])), 'same'))

def get_corr (im, pars):

    # transform 2-d image to 1-d total signal vectors
    totsig, totsigla = get_totsig(im,pars.la)

    ln1 = im.npix / 2
    ln2 = im.npix
    quad_len = ln1*ln1

    # compute staypuft signal for each pixel
    mask = totsig > 40.0
    p0 = np.fabs(pars.ampscale * mask * totsig)
    p1 = np.fabs(pars.ampscale * mask * totsigla)

    ekern = np.exp(-np.arange(ln1*pars.tx)/pars.hh)
    qkern = pars.a1*np.arange(ln1*pars.tx) + pars.a2

    e = convolve (p0, ekern, mode='full')
    q = convolve (p1, qkern, mode='full')
    b = e[0:quad_len] + q[0:quad_len]

    # transform the correction vector back into a 2-d image quad
    b = b[::-1]
    b = np.reshape (b, (ln1,ln1))
    b = np.transpose(b)

    # replicate the correction into all 4 full image quads
    im.data[0:ln1,0:ln1] = b
    im.data[0:ln1,ln1:ln2] = b
    im.data[ln1:ln2,0:ln1] = b
    im.data[ln1:ln2,ln1:ln2] = b

    return im

 def test_consistency_convolve_funcs(self):
     # Compare np.convolve, signal.convolve, signal.convolve2d
     a = np.arange(5)
     b = np.array([3.2, 1.4, 3])
     for mode in ["full", "valid", "same"]:
         assert_almost_equal(np.convolve(a, b, mode=mode), signal.convolve(a, b, mode=mode))
         assert_almost_equal(np.squeeze(signal.convolve2d([a], [b], mode=mode)), signal.convolve(a, b, mode=mode))

def basefreq(audiofile):
    """
    This function reads in the audio file and does the hann windowed fft of 
    the right input. It then smooths the output using a gaussian filter and
    then finds the peaks. It returns the peaks in the right audio channel since
    testing showed there was no significant difference in the two.
    """
    #read the data into an ndarray using scikits-audiolab        
    data, rate, enc = al.aiffread(audiofile)
    #split the left and right channel
    datar = data[:,1]
    datal = data[:,0]
    #take the fft of both of the channels with the hann window applied
    #the hann window reduces spectral leakage in the FFT     
    dftr = abs(fft.fft(datar*signal.hann(len(datar))))
    dftl = abs(fft.fft(datal*signal.hann(len(datal))))
    #compute the frequencies in the FFT
    freq = float(rate)/float(len(datar))
    freqs = np.arange(len(dftr)/2+99)*freq
    dftr = dftr[0:np.size(dftr)/2]
    dftl = dftl[0:np.size(dftr)/2]
    #smooth the fft with a gaussian
    c = signal.gaussian(100,20)
    dftr = signal.convolve(dftr,c)
    dftl = signal.convolve(dftl,c)
    #find the significant peaks in each channel
    peaksr = findpeaks(dftr,freqs)
    peaksl = findpeaks(dftl,freqs)
    #plot the output fft for testing
    #plt.plot(freqs,dftr)
    #plt.show()
    #print peaksr
    return peaksr

def compute_harris_response(image):
    """ compute the Harris corner detector response function 
        for each pixel in the image"""

    #derivatives
    imx, imy = filtertools.gauss_derivatives(image, 3)

    #kernel for blurring
    gauss = filtertools.gauss_kernel(3)

    #compute components of the structure tensor
    Wxx = signal.convolve(imx*imx, gauss, mode='same')
    Wxy = signal.convolve(imx*imy, gauss, mode='same')
    Wyy = signal.convolve(imy*imy, gauss, mode='same')

    #determinant and trace
    Wdet = Wxx*Wyy - Wxy**2
    Wtr = Wxx + Wyy

    if numpy.count_nonzero(Wtr) == 0:
        return
    print Wtr.shape
    print numpy.count_nonzero(Wtr)

    return Wdet / Wtr

def ddwt(x, num_scales):
    """
    Дискретное вейвлет-преобразование без прореживания
    :param x: входной сигнал
    :param num_scales: число уровней разложения
    :return:
    """
    h = np.array([1, 3, 3, 1], float) / 8
    g = np.array([2, -2], float)
    signal_len = len(x)

    detail = []
    approx = []
    ap = x.copy()
    detail.append(ap.copy())  # на нулевом уровне храним исходный сигнал
    approx.append([])

    for s in range(num_scales):
        dly = 2**s
        hif = convolve(ap, g, mode="full")[dly:dly+signal_len]
        detail.append(hif)
        approx.append(ap)
        if s < num_scales-1:
            ap = convolve(ap, h, mode="full")[dly:dly+signal_len]
            dly_lo = len(h)-1
            ap[:dly_lo] = ap[dly_lo] # хак
            # вместо прореживания сигналов (Маллат) расширяем характеристики
            # фильтров
            h = fexpand(h)
            g = fexpand(g)

    return approx, detail

def lsf_to_lpc(all_lsf):
    if len(all_lsf.shape) < 2:
        all_lsf = all_lsf[None]
    order = all_lsf.shape[1]
    all_lpc = np.zeros((len(all_lsf), order + 1))
    for i in range(len(all_lsf)):
        lsf = all_lsf[i]
        zeros = np.exp(1j * lsf)
        sum_zeros = zeros[::2]
        diff_zeros = zeros[1::2]
        sum_zeros = np.hstack((sum_zeros, np.conj(sum_zeros)))
        diff_zeros = np.hstack((diff_zeros, np.conj(diff_zeros)))
        sum_filt = np.poly(sum_zeros)
        diff_filt = np.poly(diff_zeros)

        if order % 2 != 0:
            deconv_diff = sg.convolve(diff_filt, [1, 0, -1])
            deconv_sum = sum_filt
        else:
            deconv_diff = sg.convolve(diff_filt, [1, -1])
            deconv_sum = sg.convolve(sum_filt, [1, 1])

        lpc = .5 * (deconv_sum + deconv_diff)
        # Last coefficient is 0 and not returned
        all_lpc[i] = lpc[:-1]
    return np.squeeze(all_lpc)

    def convolve(self, f):
        from scipy.signal import convolve

        e_field = self.E_field_as_numpy()
        f_data = np.zeros((self.dim_x(), self.dim_y()))
        result = np.zeros((self.numberEnergies(), self.dim_x(), self.dim_y(), 2), dtype=np.complex128)

        print("Convolve I_X ", end="")
        for i_x, x_cooridinate in enumerate(self.absolute_x_coordinates()):
            if i_x%100 ==0:
                print(" ",i_x , end="")
            for i_y, y_cooridinate in enumerate(self.absolute_y_coordinates()):
                f_data[i_x, i_y] = f(x_cooridinate,y_cooridinate)

        for index_energy in range(self.numberEnergies()):
            for pol in (0,1):
                print("Convolving pol", pol)
                #r = convolve(f_data, f_data,mode='same')
                r = convolve(e_field[index_energy,:,:,pol].real, f_data,mode='same')
                r = r + 1j *convolve(e_field[index_energy,:,:,pol].imag, f_data,mode='same')

                for i_x, x_cooridinate in enumerate(self.absolute_x_coordinates()):
                    for i_y, y_cooridinate in enumerate(self.absolute_y_coordinates()):
                        result[index_energy, i_x , i_y , pol] = r[i_x,i_y]

        convolved_wavefront = NumpyWavefront(result, self.x_start(), self.x_end(), self.y_start(), self.y_end(), self.z())

        return convolved_wavefront

    def __init__(self, shapein, kernel, mode="full", fft=False, **kwargs):
        if fft:
            from scipy.signal import fftconvolve as convolve
            # check kernel shape parity
            if np.any(np.asarray(kernel.shape) % 2 != 1):
                raise ValueError("Kernels with non-even shapes are not handled for now.")
        else:
            from scipy.signal import convolve

        self.kernel = kernel
        self.mode = mode

        # reverse kernel
        s = (slice(None, None, -1), ) * kernel.ndim
        self.rkernel = kernel[s]
        # reverse mode
        if mode == 'full':
            self.rmode = 'valid'
        elif mode == 'valid':
            self.rmode = 'full'
        elif mode == 'same':
            self.rmode = 'same'
        # shapeout
        if mode == 'full':
            shapeout = [s + ks - 1 for s, ks in zip(shapein, kernel.shape)]
        if mode == 'valid':
            shapeout = [s - ks + 1 for s, ks in zip(shapein, kernel.shape)]
        if mode == 'same':
            shapeout = shapein

        matvec = lambda x: convolve(x, self.kernel, mode=self.mode)
        rmatvec = lambda x: convolve(x, self.rkernel, mode=self.rmode)
        NDOperator.__init__(self, shapein, shapeout, matvec, rmatvec, **kwargs)

    def test_input_swapping(self):
        small = arange(8).reshape(2, 2, 2)
        big = 1j * arange(27).reshape(3, 3, 3)
        big += arange(27)[::-1].reshape(3, 3, 3)

        out_array = array(
            [[[0 + 0j, 26 + 0j, 25 + 1j, 24 + 2j],
              [52 + 0j, 151 + 5j, 145 + 11j, 93 + 11j],
              [46 + 6j, 133 + 23j, 127 + 29j, 81 + 23j],
              [40 + 12j, 98 + 32j, 93 + 37j, 54 + 24j]],

             [[104 + 0j, 247 + 13j, 237 + 23j, 135 + 21j],
              [282 + 30j, 632 + 96j, 604 + 124j, 330 + 86j],
              [246 + 66j, 548 + 180j, 520 + 208j, 282 + 134j],
              [142 + 66j, 307 + 161j, 289 + 179j, 153 + 107j]],

             [[68 + 36j, 157 + 103j, 147 + 113j, 81 + 75j],
              [174 + 138j, 380 + 348j, 352 + 376j, 186 + 230j],
              [138 + 174j, 296 + 432j, 268 + 460j, 138 + 278j],
              [70 + 138j, 145 + 323j, 127 + 341j, 63 + 197j]],

             [[32 + 72j, 68 + 166j, 59 + 175j, 30 + 100j],
              [68 + 192j, 139 + 433j, 117 + 455j, 57 + 255j],
              [38 + 222j, 73 + 499j, 51 + 521j, 21 + 291j],
              [12 + 144j, 20 + 318j, 7 + 331j, 0 + 182j]]])

        assert_array_equal(convolve(small, big, 'full'), out_array)
        assert_array_equal(convolve(big, small, 'full'), out_array)
        assert_array_equal(convolve(small, big, 'same'),
                           out_array[1:3, 1:3, 1:3])
        assert_array_equal(convolve(big, small, 'same'),
                           out_array[0:3, 0:3, 0:3])
        assert_raises(ValueError, convolve, small, big, 'valid')
        assert_array_equal(convolve(big, small, 'valid'),
                           out_array[1:3, 1:3, 1:3])

def filtervertical(H, minleafDomain):
	"""
	This function applies the matched filter on the horizontal 2D histogram.
	Returns filter response for both directions, as stairs can face both ways.
	https://www.youtube.com/watch?v=S7qbelm_4Y8 --> explains matched filter good, couldn't make spectralpython matched filter work
	"""
	# from skimage.feature import canny
	# edges = canny(H)
	# plt.imshow(edges,cmap=plt.cm.gray)
	# plt.show()

	filt = create_matched_filter(minleafDomain)
	# Note that the convolution of the time-reversed wavelet is identical to cross-correlation of the wavelet with the wavelet (autocorrelation) in the input signal --> http://crewes.org/ForOurSponsors/ResearchReports/2002/2002-46.pdf
	
	filty=np.transpose(filt) # transpose to also get stairs in other direction
	fr1=signal.convolve(H,filt, mode='same')

	# plt.subplot(1,3,1)
	# plt.imshow(fr1,cmap='spectral',interpolation='none')
	# plt.title('vert matched filter')
	# plt.colorbar()

	fr2=signal.convolve(H,filty, mode='same')
	# fr3=signal.convolve2d(H,filty, mode='same') # should givve the same result as fr2
	# fr3 = fr3**2 # doesn't really work as there are very high 'outliers' that just take everything away

	# plt.subplot(1,3,2)
	# plt.imshow(fr2,cmap='spectral',interpolation='none')
	# plt.title('vert matched filter transpose')
	# plt.colorbar()

	return fr1, fr2

def gaussian( u, v, size) :
    """Smooths the velocity field with a Gaussian kernel.
    
    Parameters
    ----------
    u : 2d np.ndarray
        the u velocity component field
        
    v : 2d np.ndarray
        the v velocity component field
        
    size : int
        the half width of the kernel. Kernel
        has shape 2*size+1
        
    Returns
    -------
    uf : 2d np.ndarray
        the smoothed u velocity component field
        
    vf : 2d np.ndarray
        the smoothed v velocity component field    
        
    """
    g = _gaussian_kernel( size=size )
    uf = convolve( u, g, mode='same')
    vf = convolve( v, g, mode='same')
    return uf, vf

def computeTensor(im, sigmaG=1, factorSigma=4):

    # returns 3d array of the size of the image
    # yx stores xx, yx, and yy components of the tensor

    # get the luminance of the image, use [0.3, 0.6, 0.1]
    # use numpy's dot
    # blur the image
    imLum = lum(im)
    
    imLumBlurred = zeros( ( height(im), width(im) ) )

    ndimage.filters.gaussian_filter( imLum, sigmaG, 0, imLumBlurred )

    gradX = signal.convolve(imLumBlurred, Sobel, mode='same')
    gradY = signal.convolve(imLumBlurred, transpose(Sobel), mode='same')
    
    # construct 3 2d arrays of the elements of the tensor
    gradXX = gradX*gradX
    gradYY = gradY*gradY
    gradXY = gradX*gradY

    ndimage.filters.gaussian_filter( gradXX, sigmaG * factorSigma, 0, gradXX )
    ndimage.filters.gaussian_filter( gradXY, sigmaG * factorSigma, 0, gradXY )
    ndimage.filters.gaussian_filter( gradYY, sigmaG * factorSigma, 0, gradYY )


    # construct RGB image based on these vals
    out = constantIm(height(im), width(im), 0.0)

    out[:,:,0] = gradXX
    out[:,:,1] = gradXY
    out[:,:,2] = gradYY

    return out

def gaussian_differentiation_kernel(sigma, num_stds, order, delta, scale):
    """
    http://en.wikipedia.org/wiki/Scale_space#Gaussian_derivatives
    :param sigma:
    :param num_stds:
    :param order:
    :param delta:
    :param scale:
    :return:
    """
    delta = list(delta)
    g = gaussian_kernel(sigma, num_stds)
    d = np.ones((1,) * len(order))
    for i in range(len(order)):
        _d = differentiation_kernel(order[i]).astype(float)
        if len(_d) % 2 != 1:
            _d = (np.pad(_d, ((0, 1),), 'constant') + np.pad(_d, ((1, 0),), 'constant')) / 2.
            delta[i] *= 2
        _d /= delta[i] ** order[i]
        _d *= np.sqrt(scale[i]) ** order[i]
        shp = np.ones(len(order), int)
        shp[i] = _d.shape[0]
        _d.shape = shp
        d = spsig.convolve(d, _d)
    d = spsig.convolve(g, d)
    return d

def find_peaks(fft_set, sign="+", alpha=0.15, threshold=2.0):
    """
    Method to find peaks on a fft set.
    
    :type threshold: object
    :param fft_set:
    :param sign:
    :param alpha: noise threshold
    :param threshold:
    :return:
    """

    if sign == "-":
        fft_set = -fft_set

    # Get derivative
    derivation_vector = [1, 0, -1]
    d_fft_set = convolve(fft_set, derivation_vector, "same")

    # Checking for sign-flipping and derivative
    _sign = np.sign(d_fft_set)
    d_sign = convolve(_sign, derivation_vector, "valid")

    candidates = np.where(d_fft_set > 0)[0] + (len(derivation_vector) - 1)

    peaks = sorted(set(candidates).intersection(np.where(d_sign == -2)[0] + 1))

    # Noise remover
    peaks = np.array(peaks)[fft_set[peaks] > alpha]

    return clean_adjacent_points(peaks, fft_set, float(threshold))

def makePriorGrid(zs,rs,dr,outfile,rmin=15,rmax=30):
    rMids = np.arange(16.,28.,0.1)
    drs = 0.5*np.ones(len(rMids))
    drs[rMids<23]=0.5
    drs[rMids<20]=1.0
    drs[rMids<19]=2.0
    zEdges=np.arange(-0.075,5.075,0.1)
    zMids = (zEdges[1:]+zEdges[:-1])/2.
    allH = []
    for i in range(len(rMids)):
        #print(rMids[i])
        rMsk = np.ma.masked_outside(rs,rMids[i]-dr,rMids[i]+dr)
        zMsk = np.ma.masked_array(zs,mask=np.ma.getmask(rMsk)).compressed()

        h = np.histogram(zMsk,bins=zEdges)[0]
        kernel=np.ones(5)*(1./5.)
        h2=sig.convolve(h,kernel,mode='same')
        h3=sig.convolve(h2,kernel,mode='same')
        g = interp1d(zMids,h3,bounds_error=False, fill_value=0.0)
        tot = integrate.quad(g,0.,7.)
        h3 = h3/tot[0]
        
        if i%5==0:
            plt.plot(zMids,h3,lw=3,alpha=0.75,color=cm.jet(i/len(rMids)),label='r='+str(rMids[i]))
        else:
            plt.plot(zMids,h3,lw=3,alpha=0.5,color=cm.jet(i/len(rMids)))
        allH.append(h3)
    return([rMids,zMids,np.array(allH)])

  def __init__ (self, var, saxis, kernel, fft):
  # {{{
    ''' __init__()'''
    import numpy as np

    assert len(kernel) <= var.shape[saxis], 'Kernel must not be longer than dimension being smoothed.'

    # Construct new variable
    self.saxis = saxis
    self.var = var
    self.kernel = kernel
    self.fft = fft
    self.klen = len(kernel)

    # Normalize and reshape kernel
    self.kernel /= np.sum(self.kernel)
    self.kernel.shape = [self.klen if i == saxis else 1 for i in range(var.naxes)]

    # Determine which convolution function to use
    from scipy import signal as sg
    tdata = np.ones(len(kernel), 'd')
    if self.fft:
      try:
        sg.fftconvolve(tdata, kernel, 'same', old_behaviour=False)
        self._convolve = lambda x, y, z: sg.fftconvolve(x, y, z, old_behaviour=False)
      except TypeError:
        self._convolve = sg.fftconvolve
    else:
      try:
        sg.convolve(tdata, kernel, 'same', old_behaviour=False)
        self._convolve = lambda x, y, z: sg.convolve(x, y, z, old_behaviour=False)
      except TypeError:
        self._convolve = sg.convolve

    Var.__init__(self, var.axes, var.dtype, name=var.name, atts=var.atts, plotatts=var.plotatts)

def deriv2D(data,axis=-1,dx=1.0,noise_suppression=True):
  """ Takes 1D or 2D Derivative of 2D array using convolution
	
	result = deriv2D(data)
	result = deriv2D(data, dx)
	
	output is 2D (if only one axis specified)
	output is 3D if no axis specified [nx,ny,2] with the third dimension being [dfdx, dfdy]
	
	keywords:
	axis = 0/1  If no axis specified 2D derivative will be returned
	dx = 1.0    axis spacing, must be 2D if 2D deriv is taken - default is [1.0,1.0]
	noise_suppression = True   noise suppressing coefficients used to take derivative - default = True
  """
  s = data.shape
  if axis > len(s)-1:
    raise RuntimeError("ERROR: axis out of bounds for derivative")
 
  if noise_suppression:  
    if s[axis] < 11:
      raise RuntimeError("Data too small to use 11th order method")
    tmp = array([old_div(-1.0,512.0),old_div(-8.0,512.0),old_div(-27.0,512.0),old_div(-48.0,512.0),old_div(-42.0,512.0),0.0,old_div(42.0,512.0),old_div(48.0,512.0),old_div(27.0,512.0),old_div(8.0,512.0),old_div(1.0,512.0)])
  else:
    if s[axis] < 9:
      raise RuntimeError("Data too small to use 9th order method")
    tmp = array([old_div(1.0,280.0),old_div(-4.0,105.0),old_div(1.0,5.0),old_div(-4.0,5.0),0.0,old_div(4.0,5.0),old_div(-1.0,5.0),old_div(4.0,105.0),old_div(-1.0,280.0)])       
  
  N = old_div((tmp.size-1),2)
  if axis==1:
    W = transpose(tmp[:,None])
    data_deriv = convolve(data,W,mode='same')/dx*-1.0
    for i in range(s[0]):
      data_deriv[i,0:N-1] = old_div(deriv(data[i,0:N-1]),dx)
      data_deriv[i,s[1]-N:] = old_div(deriv(data[i,s[1]-N:]),dx)
    
  elif axis==0:
    W = tmp[:,None]
    data_deriv = convolve(data,W,mode='same')/dx*-1.0
    for i in range(s[1]):
      data_deriv[0:N-1,i] = old_div(deriv(data[0:N-1,i]),dx)
      data_deriv[s[0]-N:,i] = old_div(deriv(data[s[0]-N:,i]),dx)
  else:
    data_deriv = zeros((s[0],s[1],2))
    if len(dx)==1:
      dx = array([dx,dx])
    
    W = tmp[:,None]#transpose(multiply(tmp,ones((s[1],tmp.size))))
    data_deriv[:,:,0] = convolve(data,W,mode='same')/dx[0]*-1.0
    for i in range(s[1]):
      data_deriv[0:N-1,i,0]  =  old_div(deriv(data[0:N-1,i]),dx[0])
      data_deriv[s[0]-N:s[0]+1,i,0] = old_div(deriv(data[s[0]-N:s[0]+1,i]),dx[0])
    
    W = transpose(tmp[:,None])#multiply(tmp,ones((s[0],tmp.size)))
    data_deriv[:,:,1] = convolve(data,W,mode='same')/dx[1]*-1.0
    for i in range(s[0]):
      data_deriv[i,0:N-1,1] = old_div(deriv(data[i,0:N-1]),dx[1])
      data_deriv[i,s[1]-N:s[1]+1,1] = old_div(deriv(data[i,s[1]-N:s[1]+1]),dx[1])

  return data_deriv

def kernel_func(X, Y):
    r"""Smooth histogram Y with kernel Y

    """
    if Y is None:
        return X

    return (convolve(X, Y, mode="same")
            / convolve(np.ones_like(X), Y, mode="same"))