def test_default_level():
    # default level is the maximum permissible for the transformed axes
    data = np.ones((128, 32, 4))
    wavelet = ('db8', 'db1')
    for dec_func in [pywt.wavedec2, pywt.wavedecn]:
        for axes in [(0, 1), (2, 1), (0, 2)]:
            c = dec_func(data, wavelet, axes=axes)
            max_lev = np.min([pywt.dwt_max_level(data.shape[ax], wav)
                              for ax, wav in zip(axes, wavelet)])
            assert_equal(len(c[1:]), max_lev)

    for ax in [0, 1]:
        c = pywt.wavedecn(data, wavelet[ax], axes=(ax, ))
        assert_equal(len(c[1:]),
                     pywt.dwt_max_level(data.shape[ax], wavelet[ax]))

def my_waverec(data, wavelet, mode='sym', level=None):
    """
    Multilevel 1D Discrete Wavelet Transform of data.
    Returns coefficients list - [cAn, cDn, cDn-1, ..., cD2, cD1]

    data    - input data
    wavelet - wavelet to use (Wavelet object or name string)
    mode    - signal extension mode, see MODES
    level   - decomposition level. If level is None then it will be
              calculated using `dwt_max_level` function.
    """

    if not isinstance(wavelet, pywt.Wavelet):
        wavelet = pywt.Wavelet(wavelet)

    if level is None:
        level = pywt.dwt_max_level(len(data), wavelet.dec_len)
    elif level < 0:
        raise ValueError("Level value of %d is too low . Minimum level is 0." % level)

    a = data[0:2]
    d = data[2:4]
    for i in xrange(int(level)):
        offs = pow(2,i+2)
        a = pywt.idwt(a, d, wavelet, mode)
        d = data[offs:2*offs]

    return a

def test_wavedecn_coeff_reshape_even():
    # verify round trip is correct:
    #   wavedecn - >coeffs_to_array-> array_to_coeffs -> waverecn
    # This is done for wavedec{1, 2, n}
    rng = np.random.RandomState(1234)
    params = {'wavedec': {'d': 1, 'dec': pywt.wavedec, 'rec': pywt.waverec},
              'wavedec2': {'d': 2, 'dec': pywt.wavedec2, 'rec': pywt.waverec2},
              'wavedecn': {'d': 3, 'dec': pywt.wavedecn, 'rec': pywt.waverecn}}
    N = 28
    for f in params:
        x1 = rng.randn(*([N] * params[f]['d']))
        for mode in pywt.Modes.modes:
            for wave in wavelist:
                w = pywt.Wavelet(wave)
                maxlevel = pywt.dwt_max_level(np.min(x1.shape), w.dec_len)
                if maxlevel == 0:
                    continue

                coeffs = params[f]['dec'](x1, w, mode=mode)
                coeff_arr, coeff_slices = pywt.coeffs_to_array(coeffs)
                coeffs2 = pywt.array_to_coeffs(coeff_arr, coeff_slices,
                                               output_format=f)
                x1r = params[f]['rec'](coeffs2, w, mode=mode)

                assert_allclose(x1, x1r, rtol=1e-4, atol=1e-4)

def my_wavedec(data, wavelet, mode='sym', level=None):
    """
    Multilevel 1D Discrete Wavelet Transform of data.
    Returns coefficients list - [cAn, cDn, cDn-1, ..., cD2, cD1]

    data    - input data
    wavelet - wavelet to use (Wavelet object or name string)
    mode    - signal extension mode, see MODES
    level   - decomposition level. If level is None then it will be
              calculated using `dwt_max_level` function.
    """

    if not isinstance(wavelet, pywt.Wavelet):
        wavelet = pywt.Wavelet(wavelet)

    if level is None:
        level = pywt.dwt_max_level(len(data), wavelet.dec_len)
    elif level < 0:
        raise ValueError("Level value of %d is too low . Minimum level is 0." % level)

    coeffs_list = []

    a = data
    for i in xrange(level):
        a, d = pywt.dwt(a, wavelet, mode)
        d = list(d)
        #d.reverse()
        coeffs_list.append(d)

    a = list(a)
    #a.reverse()
    coeffs_list.append(a)
    #coeffs_list.reverse()

    return coeffs_list

	def filterData(self, icurr, Fs):
		"""
			Denoise an ionic current time-series and store it in self.eventData

			:Parameters:
				- `icurr` :	ionic current in pA
				- `Fs` :	original sampling frequency in Hz
		"""
		# self.eventData=icurr
		self.Fs=Fs

		# Set up the wavelet
		w=pywt.Wavelet(self.waveletType)

		# Calculate the maximum wavelet level for the data length
		self.maxWaveletLevel=pywt.dwt_max_level(len(icurr), filter_len=w.dec_len)

		# Perform a wavelet decomposition to the specified level
		wcoeff = pywt.wavedec(icurr, w, mode='sym', level=self.waveletLevel)

		# Perform a simple threshold by setting all the detailed coefficients
		# up to level n-1 to zero
		thresh=np.std(wcoeff[-1])*self._thselect(wcoeff, self.waveletThresholdSubType)
		thrfunc=self.thrtypedict[self.waveletThresholdType]

		# print thresh, np.std(wcoeff[-1])
		wcoeff[1:] = [ thrfunc(wc, thresh) for wc in wcoeff[1:] ]
		# for i in range(1, self.waveletLevel):
		# 	wcoeff[-i]=np.zeros(len(wcoeff[-i]))

		# Reconstruct the signal with the thresholded wavelet coefficients
		self.eventData = pywt.waverec(wcoeff, self.waveletType, mode='sym')

def c_dists(Y,use_swt=True,level_weights=False):
	w = pywt.Wavelet('sym2')
	if use_swt:
		L = pywt.swt_max_level(Y.shape[0])
		C = [pywt.swt(Y[:,i],w,level=L) for i in range(Y.shape[1])]
		C = [[list(reshape(l[0],-1)) + list(reshape(l[1],-1)) for l in c] for c in C]
	else:
		L = pywt.dwt_max_level(Y.shape[0],w)
		C = [pywt.wavedec(Y[:,i],w,level=L) for i in range(Y.shape[1])]
	if level_weights:
		if use_swt:
			raise NameError('No level weights with SWT')
		Wc = [1. for x in range(1,L+1)]
		D = zeros((len(C),len(C)))
		for i in range(len(C)):
			for j in range(i+1,len(C)):
				d = sum([distance.cosine(C[i][x],C[j][x])*Wc[x] for x in range(L)])/sum(Wc)
				D[i,j] = d
				D[j,i] = d
		return D
	else:
		Cn = []
		for c in C:
			cn = []
			for l in c:
				cn += list(l)
			Cn.append(cn)
		return abs(pdist(Cn,'cosine'))

def wp(S, costf, wavelet="db4", mode=pywt.MODES.ppd, level=None):
    '''
    Returns the 1D discrete wavelet packet transformation, with the best basis according
    to the given cost function, for the given 1D input signal.
    @param S:         Input signal.
                      Both single and double precision floating-point data types are supported
                      and the output type depends on the input type. If the input data is not
                      in one of these types it will be converted to the default double precision
                      data format before performing computations.
    @param costf:      The (single parameter) cost function that must be used while
                      searching for the best basis.
    @param wavelet:   Wavelet to use in the transform. 
                      This must be a name of the wavelet from the wavelist() list.
    @param mode:      Signal extension mode to deal with the border distortion problem.
                      The default mode is periodic-padding.
    @param level:     Number of decomposition steps to perform. If the level is None, then the
                      full decomposition up to the level computed with dwt_max_level() function for
                      the given data and wavelet lengths is performed.
    @return:          A list containing the nodes of the 1D discrete wavelet packet transformation,
                      with the best basis according to the given cost function, for the given input signal. 
    '''
    if (level == None):
        level = pywt.dwt_max_level(S.shape[0], pywt.Wavelet(wavelet))
    
    #Data collection step
    Nodes = collect(S, wavelet=wavelet, mode=mode, level=level)
    #Dynamic programming upstream traversal
    mark(Nodes, costf)
    #node.print_nodes(Nodes)
    #Dynamic programming downstream traversal
    Result = []
    traverse(Nodes[0][0], Nodes, Result)
    traverse(Nodes[0][1], Nodes, Result)
    return sorted(Result, cmp=node.compare_low_level_first, reverse=False)

def test_wavelet_denoising_levels():
    rstate = np.random.RandomState(1234)
    ndim = 2
    N = 256
    wavelet = 'db1'
    # Generate a very simple test image
    img = 0.2*np.ones((N, )*ndim)
    img[[slice(5, 13), ] * ndim] = 0.8

    sigma = 0.1
    noisy = img + sigma * rstate.randn(*(img.shape))
    noisy = np.clip(noisy, 0, 1)

    denoised = restoration.denoise_wavelet(noisy, wavelet=wavelet)
    denoised_1 = restoration.denoise_wavelet(noisy, wavelet=wavelet,
                                             wavelet_levels=1)
    psnr_noisy = compare_psnr(img, noisy)
    psnr_denoised = compare_psnr(img, denoised)
    psnr_denoised_1 = compare_psnr(img, denoised_1)

    # multi-level case should outperform single level case
    assert_(psnr_denoised > psnr_denoised_1 > psnr_noisy)

    # invalid number of wavelet levels results in a ValueError
    max_level = pywt.dwt_max_level(np.min(img.shape),
                                   pywt.Wavelet(wavelet).dec_len)
    assert_raises(ValueError, restoration.denoise_wavelet, noisy,
                  wavelet=wavelet, wavelet_levels=max_level+1)
    assert_raises(ValueError, restoration.denoise_wavelet, noisy,
                  wavelet=wavelet, wavelet_levels=-1)

def wavelet_trans(l, keep_proportion = 0.5):
    #TODO binning per decomposition level?
    #TODO explain why db4
    if not (l is None or l == []):
        coeffs = pywt.wavedec(l, 'db4', level=pywt.dwt_max_level(len(l),pywt.Wavelet('db4')))
        return merge(*coeffs[0:int(len(coeffs)*keep_proportion)])
    else:
        return l

    def full_coeff_len(datalen, filtlen, mode):
        max_level = wt.dwt_max_level(datalen, filtlen)
        total_len = 0

        for i in xrange(max_level):
            datalen = wt.dwt_coeff_len(datalen, filtlen, mode)
            total_len += datalen

        return total_len + datalen

def test_wavedec():
    x = [3, 7, 1, 1, -2, 5, 4, 6]
    db1 = pywt.Wavelet('db1')
    cA3, cD3, cD2, cD1 = pywt.wavedec(x, db1)
    assert_almost_equal(cA3, [8.83883476])
    assert_almost_equal(cD3, [-0.35355339])
    assert_allclose(cD2, [4., -3.5])
    assert_allclose(cD1, [-2.82842712, 0, -4.94974747, -1.41421356])
    assert_(pywt.dwt_max_level(len(x), db1) == 3)

def startup():
    global w, N, L0, mode, level

    print pywt.families()
    print pywt.wavelist('db')

    w = pywt.Wavelet('db6')
    mode = pywt.MODES.per

    print w
    print "vanishing_moments_psi:", w.vanishing_moments_psi
    print "vanishing_moments_phi:", w.vanishing_moments_phi

    N = 2**9
    print "max level = ", pywt.dwt_max_level(N, w.dec_len)
    L0 = numpy.zeros((N,N), 'double')
    if True:
        for i in xrange(0,N):
            L0[i][i-1], L0[i][i], L0[i-1][i] = (1., -2., 1.)
    else:
        for i in xrange(1,N):
            L0[i][i-1], L0[i][i], L0[i-1][i] = (1., -2., 1.)
        L0[0][0] = -2.

#L0[0][N-1], L0[0][0], L0[N-1][0] = (1, -2, 1)

#L0 = numpy.eye(N)

#for i in xrange(0,N):
        #L0[i] = [1,2,3,4,5,6,7,8]

    numpy.core.arrayprint.set_printoptions(threshold=N*N+1, linewidth=100000)

#print L0

#coeffs = pywt.wavedec2(L0, w) #, level=pywt.dwt_max_level(N, w.dec_len))
#print coeffs
#print pywt.waverec2(coeffs, w)
#
#coeffs = pywt.wavedec2(L0, w) #, level=pywt.dwt_max_level(N, w.dec_len))
#print coeffs

    print "max level = ", pywt.dwt_max_level(N, w.dec_len)

def wavelet_levels(Y):
	w = pywt.Wavelet('sym2')
	levels = pywt.dwt_max_level(Y.shape[0],w)
	w0 = pywt.wavedec(Y[:,0],w,level=levels)[1:]
	L = [np.empty((Y.shape[1],len(x))) for x in w0]
	for i in range(Y.shape[1]):
		wd = pywt.wavedec(Y[:,i],w)[1:]
		for j,x in enumerate(wd):
			L[j][i,:] = x
	return L,[Y.shape[0]/len(x) for x in w0]

def haar_decomp(img_file_name):
    img = Image.open(img_file_name).convert('L')
    x,y = img.size
    coeffs = pywt.wavedec2(img, 'haar', level=pywt.dwt_max_level(max(x,y), pywt.Wavelet('haar')))
    params = coeffs[0]
    final_arr = []
    for i in xrange(1,len(coeffs)):
	    for j in coeffs[i]:
		    final_arr = np.concatenate((final_arr, np.ndarray.flatten(j)))
    return params, coeffs[1:], final_arr

def wavelet_avg(Y,X=None,reshape=False,plot_xy=False,zero_thresh=1*10**-7,Names=None,Title=None):
	if not X:
		X = range(Y.shape[0])
	w = pywt.Wavelet('sym2')
	L = pywt.dwt_max_level(Y.shape[0],w)
	#L = pywt.swt_max_level(Y.shape[0])
	Zavg = zeros((L,len(X)))
	Zavgabs = zeros((L,len(X)))
	if reshape:
		e = [0,L,0,L]
	else:
		e = [X[0],X[-1],0,L]
	sp = 1
	if not Names:
		Names = ["Series "+str(y+1) for y in range(Y.shape[1])]
	if Title:
		pylab.suptitle(Title)
	if plot_xy:
		tp = str(Y.shape[1]+3)
		pylab.subplot(int(''.join([tp,'1',str(sp)])))
		P = pylab.plot(X,Y)
		pylab.legend(P,Names,prop={'size':6})
		pylab.xlim([0,Y.shape[0]-1])
		pylab.title("Relative Abundance Time Series")
		sp += 1
	else:
		tp = str(Y.shape[1]+2)
	for y in range(Y.shape[1]):
		C = pywt.wavedec(Y[:,y],w,level=L)
		#C = pywt.swt(Y[:,y],w,level=L)
		#C = ['dummy'] + [c[1] for c in C]
		Z = wavelet_matrix(C,L,len(X))
		pylab.subplot(int(''.join([tp,'1',str(sp)])))
		pylab.imshow(abs(Z),extent=e)
		pylab.title(Names[y])
		sp += 1
		Zavg += Z
		Zavgabs += abs(Z)
	Zavg /= Y.shape[1]
	Zavg[Zavg<zero_thresh] = 0
	pylab.subplot(int(''.join([tp,'1',str(sp)])))
	pylab.imshow(abs(Zavg),extent=e)
	pylab.title("Decomposition Avg")
	sp += 1
	Zavgabs /= Y.shape[1]
	Zavgabs[Zavgabs<zero_thresh] = 0
	pylab.subplot(int(''.join([tp,'1',str(sp)])))
	pylab.imshow(abs(Zavgabs),extent=e)
	pylab.title("abs(Decomposition) Sum")
	sp += 1
	pylab.show()
	return Zavg

def test_compare_downcoef_coeffs():
    rstate = np.random.RandomState(1234)
    r = rstate.randn(16)
    # compare downcoef against wavedec outputs
    for nlevels in [1, 2, 3]:
        for wavelet in pywt.wavelist():
            wavelet = pywt.Wavelet(wavelet)
            max_level = pywt.dwt_max_level(r.size, wavelet.dec_len)
            if nlevels <= max_level:
                a = pywt.downcoef('a', r, wavelet, level=nlevels)
                d = pywt.downcoef('d', r, wavelet, level=nlevels)
                coeffs = pywt.wavedec(r, wavelet, level=nlevels)
                assert_allclose(a, coeffs[0])
                assert_allclose(d, coeffs[1])

def test_dwt_max_level():
    assert_(pywt.dwt_max_level(16, 2) == 4)
    assert_(pywt.dwt_max_level(16, 8) == 1)
    assert_(pywt.dwt_max_level(16, 9) == 1)
    assert_(pywt.dwt_max_level(16, 10) == 0)
    assert_(pywt.dwt_max_level(16, np.int8(10)) == 0)
    assert_(pywt.dwt_max_level(16, 10.) == 0)
    assert_(pywt.dwt_max_level(16, 18) == 0)

    # accepts discrete Wavelet object or string as well
    assert_(pywt.dwt_max_level(32, pywt.Wavelet('sym5')) == 1)
    assert_(pywt.dwt_max_level(32, 'sym5') == 1)

    # string input that is not a discrete wavelet
    assert_raises(ValueError, pywt.dwt_max_level, 16, 'mexh')

    # filter_len must be an integer >= 2
    assert_raises(ValueError, pywt.dwt_max_level, 16, 1)
    assert_raises(ValueError, pywt.dwt_max_level, 16, -1)
    assert_raises(ValueError, pywt.dwt_max_level, 16, 3.3)

def wavelet_filter(data, thresh, scale, dyadmin):
    ndat   = np.size(data)

    # Use the Coiflet 3 wavelet (18 data points)
    wavelet = pywt.Wavelet('coif3')

    #per - periodization - is like periodic-padding but gives the smallest possible number
    #of decomposition coefficients. IDWT must be performed with the same mode.
    mode = 'per'

    # dwt_max_level() == floor( log(data_len/(filter_len-1)) / log(2) )
    max_level = pywt.dwt_max_level(data_len=len(data), filter_len=wavelet.dec_len)

    # Perform the DWT using the the full decomposition up to the level computed
    # with dwt_max_level() function for the given data and wavelet lengths.
    coeff = pywt.wavedec(data, wavelet, mode, level=max_level)

    # Use max line width as criterion for removing wavelet
    # coefficients with features on smaller scales
    width = np.where(thresh > 2, thresh, 2)  # this is in points

    # find next power of 2
    pow2 = 0
    while width > 2**pow2: pow2 = pow2 + 1

    # figure out which indices in the wavelet array to zero
    dyadmin = int(np.log(ndat)/np.log(2)) - dyadmin            # Hard set min
    dyadlw  = int(np.log(ndat)/np.log(2)) - pow2 - scale       # LW multiplier set
    windex  = np.where(dyadlw < dyadmin, dyadlw, dyadmin)

    # based on coarsest value for the given mother wavelet, get min set of coeffs
    minwindex = int(thresh)
    windex = np.where(windex > minwindex, windex, minwindex)

    # Calculate the number of coefficients that need to be zeroed (windex)
    windex = int(np.log2(len(data))) - windex

    # Zero the last 'windex' coefficents where windex is the number of zeroed coefficients
    ncoeff = len(coeff)
    for i in range(ncoeff - windex, ncoeff):
        coeff[i] = 0.0 * np.array(coeff[i])

    res = pywt.waverec(coeff, wavelet, mode)

    return res

def test_wavelet_denoising_levels():
    rstate = np.random.RandomState(1234)
    ndim = 2
    N = 256
    wavelet = 'db1'
    # Generate a very simple test image
    img = 0.2*np.ones((N, )*ndim)
    img[(slice(5, 13), ) * ndim] = 0.8

    sigma = 0.1
    noisy = img + sigma * rstate.randn(*(img.shape))
    noisy = np.clip(noisy, 0, 1)

    with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
        denoised = restoration.denoise_wavelet(noisy, wavelet=wavelet)
    denoised_1 = restoration.denoise_wavelet(noisy, wavelet=wavelet,
                                             wavelet_levels=1)
    psnr_noisy = compare_psnr(img, noisy)
    psnr_denoised = compare_psnr(img, denoised)
    psnr_denoised_1 = compare_psnr(img, denoised_1)

    # multi-level case should outperform single level case
    assert_(psnr_denoised > psnr_denoised_1 > psnr_noisy)

    # invalid number of wavelet levels results in a ValueError or UserWarning
    max_level = pywt.dwt_max_level(np.min(img.shape),
                                   pywt.Wavelet(wavelet).dec_len)
    if Version(pywt.__version__) < '1.0.0':
        # exceeding max_level raises a ValueError in PyWavelets 0.4-0.5.2
        with testing.raises(ValueError):
            with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
                restoration.denoise_wavelet(
                    noisy, wavelet=wavelet, wavelet_levels=max_level + 1)
    else:
        # exceeding max_level raises a UserWarning in PyWavelets >= 1.0.0
        with expected_warnings([
                'all coefficients will experience boundary effects']):
            restoration.denoise_wavelet(
                noisy, wavelet=wavelet, wavelet_levels=max_level + 1)

    with testing.raises(ValueError):
        with expected_warnings([PYWAVELET_ND_INDEXING_WARNING]):
            restoration.denoise_wavelet(
                noisy,
                wavelet=wavelet, wavelet_levels=-1)

def wavedec(data, wavelet, mode='sym', level=None):

	if level is None:
		level = dwt_max_level(len(data), wavelet.dec_len)
	elif level < 0:
			raise ValueError( "Level value of %d is too low . Minimum level is 0." % level)

	coeffs_list = []

	a = data
	for i in xrange(level):
		a, d = dwt(a, wavelet, mode)
		coeffs_list.append(d)

	coeffs_list.append(a)
	coeffs_list.reverse()

	return coeffs_list