def test_swt_decomposition():
    x = [3, 7, 1, 3, -2, 6, 4, 6]
    db1 = pywt.Wavelet('db1')
    (cA3, cD3), (cA2, cD2), (cA1, cD1) = pywt.swt(x, db1, level=3)
    expected_cA1 = [7.07106781, 5.65685425, 2.82842712, 0.70710678,
                    2.82842712, 7.07106781, 7.07106781, 6.36396103]
    assert_allclose(cA1, expected_cA1)
    expected_cD1 = [-2.82842712, 4.24264069, -1.41421356, 3.53553391,
                    -5.65685425, 1.41421356, -1.41421356, 2.12132034]
    assert_allclose(cD1, expected_cD1)
    expected_cA2 = [7, 4.5, 4, 5.5, 7, 9.5, 10, 8.5]
    assert_allclose(cA2, expected_cA2, rtol=tol_double)
    expected_cD2 = [3, 3.5, 0, -4.5, -3, 0.5, 0, 0.5]
    assert_allclose(cD2, expected_cD2, rtol=tol_double, atol=1e-14)
    expected_cA3 = [9.89949494, ] * 8
    assert_allclose(cA3, expected_cA3)
    expected_cD3 = [0.00000000, -3.53553391, -4.24264069, -2.12132034,
                    0.00000000, 3.53553391, 4.24264069, 2.12132034]
    assert_allclose(cD3, expected_cD3)

    # level=1, start_level=1 decomposition should match level=2
    res = pywt.swt(cA1, db1, level=1, start_level=1)
    cA2, cD2 = res[0]
    assert_allclose(cA2, expected_cA2, rtol=tol_double)
    assert_allclose(cD2, expected_cD2, rtol=tol_double, atol=1e-14)

    coeffs = pywt.swt(x, db1)
    assert_(len(coeffs) == 3)
    assert_(pywt.swt_max_level(len(x)) == 3)

def _collect_coefficients(data, genome, loci):
    l = loci
    scale = genome[l.scale]
    (temps, loads) = (data['Temperature'], data['Load'])
    loads_coeffs = pywt.swt(loads, 'haar', level=scale)
    temps_coeffs = pywt.swt(temps, 'haar', level=scale)

    # The first 2^scale datapoints cannot be used to predict because of lack of
    # history. scale+1 because of the smooth array.
    a = np.zeros((len(loads) - 2**scale, 2*(scale+1)*genome[l.Aj]))

    # Collect coefficients for each scale + smooth array.
    for i in range(len(a)):
        row = []
        # cAn, the smoothest of the smooth arrays.
        for k in range(1, genome[l.Aj]+1):
            row.append(loads_coeffs[-1][0][2**scale + i - 2**scale*(k-1)])
            row.append(temps_coeffs[-1][0][2**scale + i - 2**scale*(k-1)])
        # cD, the details.
        for j in range(1, scale+1):
            for k in range(1, genome[l.Aj]+1):
                row.append(loads_coeffs[j-1][1][2**scale + i - 2**j*(k-1)])
                row.append(temps_coeffs[j-1][1][2**scale + i - 2**j*(k-1)])

        a[i] = np.array(row)

    return a

def muunna(sarake):
    aa = [];
    bb = []
    #pituus
    # floor (pituus/16000)
    loops = range(1,int(math.floor(len(sarake)/16000)))
    max_loops = int(len(sarake) - math.floor(len(sarake)/16000)*16000)
    
    wave2 = pywt.Wavelet('db4')
    
    sarake1 = np.array(sarake[0:16000])

    (cA6, cD6), (cA5, cD5), (cA4, cD4), (cA3, cD3), (cA2, cD2), (cA1, cD1) = pywt.swt(sarake1, wave2, level=6)
    aa = np.column_stack((cD1,cD2,cD3,cD4,cD5,cD6))
    
    for loop in loops:
    
        # loop yli len(col) mod 16000 tai 8000
        sarake1 = np.array(sarake[16000*loop:16000*(loop+1)])
    
        (cA6, cD6), (cA5, cD5), (cA4, cD4), (cA3, cD3), (cA2, cD2), (cA1, cD1) = pywt.swt(sarake1, wave2, level=6)
        bb = np.column_stack((cD1,cD2,cD3,cD4,cD5,cD6))
        aa = np.vstack((aa,bb))
    
    # lopussa valitaan viimeiset 16000
    # ja sijoitetaan kohtaan [(loops*16000):len(sarake)]
    sarake2 = np.array(sarake[(len(sarake) - 16000):len(sarake)])
    (cA6, cD6), (cA5, cD5), (cA4, cD4), (cA3, cD3), (cA2, cD2), (cA1, cD1) = pywt.swt(sarake2, wave2, level=6)
    bb = np.column_stack((cD1,cD2,cD3,cD4,cD5,cD6))[(len(cD1)-max_loops):len(cD1)]
    
    aa = np.vstack((aa,bb))
    #aa = np.concatenate((aa,bb), axis=0)
    #bb.append(np.concatenate(aa2, axis=0))
    return aa

def test_swt_axis():
    x = [3, 7, 1, 3, -2, 6, 4, 6]

    db1 = pywt.Wavelet('db1')
    (cA2, cD2), (cA1, cD1) = pywt.swt(x, db1, level=2)

    # test cases use 2D arrays based on tiling x along an axis and then
    # calling swt along the other axis.
    for order in ['C', 'F']:
        # test SWT of 2D data along default axis (-1)
        x_2d = np.asarray(x).reshape((1, -1))
        x_2d = np.concatenate((x_2d, )*5, axis=0)
        if order == 'C':
            x_2d = np.ascontiguousarray(x_2d)
        elif order == 'F':
            x_2d = np.asfortranarray(x_2d)
        (cA2_2d, cD2_2d), (cA1_2d, cD1_2d) = pywt.swt(x_2d, db1, level=2)

        for c in [cA2_2d, cD2_2d, cA1_2d, cD1_2d]:
            assert_(c.shape == x_2d.shape)
        # each row should match the 1D result
        for row in cA1_2d:
            assert_array_equal(row, cA1)
        for row in cA2_2d:
            assert_array_equal(row, cA2)
        for row in cD1_2d:
            assert_array_equal(row, cD1)
        for row in cD2_2d:
            assert_array_equal(row, cD2)

        # test SWT of 2D data along other axis (0)
        x_2d = np.asarray(x).reshape((-1, 1))
        x_2d = np.concatenate((x_2d, )*5, axis=1)
        if order == 'C':
            x_2d = np.ascontiguousarray(x_2d)
        elif order == 'F':
            x_2d = np.asfortranarray(x_2d)
        (cA2_2d, cD2_2d), (cA1_2d, cD1_2d) = pywt.swt(x_2d, db1, level=2,
                                                      axis=0)

        for c in [cA2_2d, cD2_2d, cA1_2d, cD1_2d]:
            assert_(c.shape == x_2d.shape)
        # each column should match the 1D result
        for row in cA1_2d.transpose((1, 0)):
            assert_array_equal(row, cA1)
        for row in cA2_2d.transpose((1, 0)):
            assert_array_equal(row, cA2)
        for row in cD1_2d.transpose((1, 0)):
            assert_array_equal(row, cD1)
        for row in cD2_2d.transpose((1, 0)):
            assert_array_equal(row, cD2)

    # axis too large
    assert_raises(ValueError, pywt.swt, x, db1, level=2, axis=5)

def get_waveletfeatures(data, w, use_dwt=True):
    #Show dwt or swt coefficients for given data and wavelet.
    w = pywt.Wavelet(w)
    a = data
    ca = []
    cd = []

    if use_dwt:
        for i in range(5):
            (a, d) = pywt.dwt(a, w, mode)
            ca.append(a)
            cd.append(d)
    else:
        coeffs = pywt.swt(data, w, 5)  # [(cA5, cD5), ..., (cA1, cD1)]
        for a, d in reversed(coeffs):
            ca.append(a)
            cd.append(d)
    
    wave_features = []
    for i in range(len(ca)): #ca1 - ca5
        c_Dsquares = []
        for j in range(len(cd[i])):
            c_Dsquares.append((cd[i][j])**2)
        c_Dsumsquares=sum(c_Dsquares)
        wave_features.append(c_Dsumsquares)
    return wave_features # Returns 

def stationary_hard_filter (y, sigma, tau, level=3):
    threshold= tau * sigma    
    coeffs = pywt.swt(y, 'db6', level)
    hcoeffs =[]
    for scale, x in enumerate(coeffs):
        hcoeffs.append(pywt.thresholding.hard(x, threshold))
    return iswt(hcoeffs, 'db6') 

def test_iswt_mixed_dtypes():
    # Mixed precision inputs give double precision output
    x_real = np.arange(16).astype(np.float64)
    x_complex = x_real + 1j*x_real
    wav = 'sym2'
    for dtype1, dtype2 in [(np.float64, np.float32),
                           (np.float32, np.float64),
                           (np.float16, np.float64),
                           (np.complex128, np.complex64),
                           (np.complex64, np.complex128)]:

        if dtype1 in [np.complex64, np.complex128]:
            x = x_complex
            output_dtype = np.complex128
        else:
            x = x_real
            output_dtype = np.float64

        coeffs = pywt.swt(x, wav, 2)
        # different precision for the approximation coefficients
        coeffs[0] = [coeffs[0][0].astype(dtype1),
                     coeffs[0][1].astype(dtype2)]
        y = pywt.iswt(coeffs, wav)
        assert_equal(output_dtype, y.dtype)
        assert_allclose(y, x, rtol=1e-3, atol=1e-3)

def test_swt_default_level_by_axis():
    # make sure default number of levels matches the max level along the axis
    wav = 'db2'
    x = np.ones((2**3, 2**4, 2**5))
    for axis in (0, 1, 2):
        sdec = pywt.swt(x, wav, level=None, start_level=0, axis=axis)
        assert_equal(len(sdec), pywt.swt_max_level(x.shape[axis]))

def test_swt_iswt_integration():
    # This function performs a round-trip swt/iswt transform test on
    # all available types of wavelets in PyWavelets - except the
    # 'dmey' wavelet. The latter has been excluded because it does not
    # produce very precise results. This is likely due to the fact
    # that the 'dmey' wavelet is a discrete approximation of a
    # continuous wavelet. All wavelets are tested up to 3 levels. The
    # test validates neither swt or iswt as such, but it does ensure
    # that they are each other's inverse.

    max_level = 3
    wavelets = pywt.wavelist()
    if 'dmey' in wavelets:
        # The 'dmey' wavelet seems to be a bit special - disregard it for now
        wavelets.remove('dmey')
    for current_wavelet_str in wavelets:
        current_wavelet = pywt.Wavelet(current_wavelet_str)
        input_length_power = int(np.ceil(np.log2(max(
            current_wavelet.dec_len,
            current_wavelet.rec_len))))
        input_length = 2**(input_length_power + max_level - 1)
        X = np.arange(input_length)
        coeffs = pywt.swt(X, current_wavelet, max_level)
        Y = pywt.iswt(coeffs, current_wavelet)
        assert_allclose(Y, X, rtol=1e-5, atol=1e-7)

def waveletDecomp(series):
    # c6, k6, k5, k4 ,k3, k2, k1 = pywt.wavedec(series,"db6",level=6)
    (c6,k6),(c5,k5),(c4,k4),(c3,k3),(c2,k2),(c1,k1)=pywt.swt(series,"haar",level=6)
    wave_matrix = []
    for i in range(10):
        #count = int(32/len(k6))
        #map(lambda n:[wave_matrix.append(n) for i in range(count)], k6)
        map(lambda n:wave_matrix.append(n), k6)
    for i in range(10):
        #count = int(32/len(k5))
        #map(lambda n:[wave_matrix.append(n) for i in range(count), k5)
        map(lambda n:wave_matrix.append(n), k5)
    for i in range(11):
        #count = int(32/len(k4))
        #map(lambda n:[wave_matrix.append(n) for i in range(count)], k4)
        map(lambda n:wave_matrix.append(n), k4)
    for i in range(11):
        #count = int(32/len(k3))
        #map(lambda n:[wave_matrix.append(n) for i in range(count)], k3)
        map(lambda n:wave_matrix.append(n), k3)
    for i in range(11):
        #count = int(32/len(k2))
        #map(lambda n:[wave_matrix.append(n) for i in range(count)], k2)
        map(lambda n:wave_matrix.append(n), k2)
    for i in range(11):
        #count = int(32/len(k1))
        #map(lambda n:[wave_matrix.append(n) for i in range(count)], k1)
        map(lambda n:wave_matrix.append(n), k1)
    return wave_matrix

def random():
    fig, (ax_w, ax_f) = plt.subplots(2)

    distributions = [.2, .4, .6, .8, .99]
    colors = ['b', 'r', 'g', 'm', 'c']

    for d,c in zip(distributions, colors):
        x_wavelet, y_wavelet, x_fourier, y_fourier = [],[],[],[]
        for i in range(int(1e3)):
            signal = np.random.rand(1e5) > d
            (cA2, cD2), (cA1, cD1) = pywt.swt(signal, 'haar', level=2)
            x_wavelet.append(np.mean(cA1))
            y_wavelet.append(np.mean(cA2))

            fourier = np.fft.fft(signal)
            freqs = np.fft.fftfreq(len(fourier))
            peaks = argrelextrema(np.abs(fourier), np.greater)

            x_fourier.append(freqs[peaks[0][0]])
            y_fourier.append(freqs[peaks[0][1]])
            
        ax_w.scatter(x_wavelet, y_wavelet, marker='o', c=c, label='w: rand(1e5) > {}'.format(d), linewidth=0)
        ax_f.scatter(x_fourier, y_fourier, marker='+', c=c, label='f: rand(1e5) > {}'.format(d))

    ax_w.set_title('Random signals')
    ax_w.set_xlabel('Wavelet transform')
    ax_f.set_xlabel('Fourier transform')
    #plt.legend()

    plt.show()

 def setup(self, n, wavelet):
     try:
         from pywt import iswt
     except ImportError:
         raise NotImplementedError("iswt not available")
     super(IswtTimeSuite, self).setup(n, wavelet)
     self.coeffs = pywt.swt(self.data, wavelet)

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 calculate_suitable_lvl(data, wv, r, swt=True):
    # stationary
    if swt:
        max_lvl = pywt.swt_max_level(len(data))

        lvl = 1
        ent = []
        pre_e = entropy(pywt.swt(data, wv, lvl), r)
        ent.append(pre_e)
        lvl += 1
        while True:
            new_e = entropy(pywt.swt(data, wv, lvl), r)
            ent.append(new_e)
            if lvl < max_lvl:
                lvl += 1
            else:
                break

        e_sorted = sorted(ent[:])
        median = e_sorted[len(e_sorted) / 2]
        lvl = ent.index(median) + 1
    # discrete
    else:
        lvl = 1
        data_e = entropy(data, r, True)
        max_lvl = pywt.dwt_max_level(len(data), wv)

        while True:
            new_e = entropy(pywt.dwt(data, wv, lvl)[lvl - 1], r, True)
            if new_e > data_e:
                break
            elif lvl == max_lvl:
                break
            else:
                lvl += 1

    if lvl == max_lvl:
        pass

    return lvl

def test_swt_decomposition():
    x = [3, 7, 1, 3, -2, 6, 4, 6]
    db1 = pywt.Wavelet("db1")
    (cA2, cD2), (cA1, cD1) = pywt.swt(x, db1, level=2)
    assert_allclose(
        cA1, [7.07106781, 5.65685425, 2.82842712, 0.70710678, 2.82842712, 7.07106781, 7.07106781, 6.36396103]
    )
    assert_allclose(
        cD1, [-2.82842712, 4.24264069, -1.41421356, 3.53553391, -5.65685425, 1.41421356, -1.41421356, 2.12132034]
    )
    expected_cA2 = [7, 4.5, 4, 5.5, 7, 9.5, 10, 8.5]
    assert_allclose(cA2, expected_cA2, rtol=1e-12)
    expected_cD2 = [3, 3.5, 0, -4.5, -3, 0.5, 0, 0.5]
    assert_allclose(cD2, expected_cD2, rtol=1e-12, atol=1e-14)

    # level=1, start_level=1 decomposition should match level=2
    res = pywt.swt(cA1, db1, level=1, start_level=1)
    cA2, cD2 = res[0]
    assert_allclose(cA2, expected_cA2, rtol=1e-12)
    assert_allclose(cD2, expected_cD2, rtol=1e-12, atol=1e-14)

    coeffs = pywt.swt(x, db1)
    assert_(len(coeffs) == 3)
    assert_(pywt.swt_max_level(len(x)) == 3)

def get_cdpp(total_ratio):
    length = total_ratio.shape[0]
    level=12
    cdpp_list = []
    for trail_length in [6, 12, 24]:
        signal_origin = np.zeros(length)
        for i in range(1500, 1500+trail_length):
            signal_origin[i] = 1.

        ratio = signal_extension(total_ratio)
        signal = signal_extension(signal_origin)

        wavelet = pywt.Wavelet('db6')

        swc_ratio = pywt.swt(ratio, wavelet, level)
        swc_signal = pywt.swt(signal, wavelet, level)
        x = []
        s = []
        for i in range(0, level):
            x.append(swc_ratio[level-1-i][1])
            s.append(swc_signal[level-1-i][1])
        x.append(swc_ratio[0][0])
        s.append(swc_signal[0][0])

        K = 50*trail_length
        sigma = variance(x, K)
        
        M = len(s)
        D = np.zeros(8192)
        for i in range(1, M+1):
            power = np.min([i, M-1])
            D += (2**(-power))*np.convolve(sigma[i-1]**(-1), s[i-1]**2, 'same')
        cdpp = (1e6)*np.sqrt(D)**(-1)
        rms_cdpp = math.sqrt(np.mean(cdpp**2))
        cdpp_list.append(rms_cdpp)
    return cdpp_list

def dwt_swt_plot(data, w, filename,DWT,mode='sp1',lev=4):
    w = pywt.Wavelet(w)
    a = data
    ca = []
    cd = []

    if DWT:
        for i in xrange(5):
            (a, d) = pywt.dwt(a, w, mode=mode)
            ca.append(a)
            cd.append(d)
    else:
        coeffs = pywt.swt(data, w) #level=lev # [(cA5, cD5), ..., (cA1, cD1)]
        for a, d in reversed(coeffs):
            ca.append(a)
            cd.append(d)
                
    pylab.figure()
    ax_main = pylab.subplot(len(ca) + 1, 1, 1)
    pylab.title(filename)
    ax_main.plot(data)
    pylab.xlim(0, len(data) - 1)
    
    for i, x in enumerate(ca):
        ax = pylab.subplot(len(ca) + 1, 2, 3 + i * 2)
        ax.plot(x, 'r')
        if DWT:
            pylab.xlim(0, len(x) - 1)
        else:
            pylab.xlim(w.dec_len * i, len(x) - 1 - w.dec_len * i)
            pylab.ylabel("A%d" % (i + 1))
                        
    for i, x in enumerate(cd):
        ax = pylab.subplot(len(cd) + 1, 2, 4 + i * 2)
        ax.plot(x, 'g')
        pylab.xlim(0, len(x) - 1)
        if DWT:
            pylab.ylim(min(0, 1.4 * min(x)), max(0, 1.4 * max(x)))
        else:  # SWT
            pylab.ylim(
                min(0, 2 * min(
                    x[w.dec_len * (1 + i):len(x) - w.dec_len * (1 + i)])),
                max(0, 2 * max(
                    x[w.dec_len * (1 + i):len(x) - w.dec_len * (1 + i)]))
                                                                )
            pylab.ylabel("D%d" % (i + 1))
            
    pylab.savefig(filename+'.pdf')

def test_swt_roundtrip_dtypes():
    # verify perfect reconstruction for all dtypes
    rstate = np.random.RandomState(5)
    wavelet = pywt.Wavelet('haar')
    for dt_in, dt_out in zip(dtypes_in, dtypes_out):
        # swt, iswt
        x = rstate.standard_normal((8, )).astype(dt_in)
        c = pywt.swt(x, wavelet, level=2)
        xr = pywt.iswt(c, wavelet)
        assert_allclose(x, xr, rtol=1e-6, atol=1e-7)

        # swt2, iswt2
        x = rstate.standard_normal((8, 8)).astype(dt_in)
        c = pywt.swt2(x, wavelet, level=2)
        xr = pywt.iswt2(c, wavelet)
        assert_allclose(x, xr, rtol=1e-6, atol=1e-7)

def denoise_waveform(wf_array, flatTimeSamples):
  #should already by a numpy array
  threshold_list = get_threshold_list()

  swt_output = pywt.swt(wf_array, wl, level=levels)
  # threshold the SWT coefficients
  apply_threshold(swt_output, 1., threshold_list)
  # inverse transform
  cA_thresh = iswt(swt_output, wl)
  
  wf_array = cA_thresh

  #re-baseline-subtract
  wf_array -= np.mean(wf_array[:flatTimeSamples])

  return wf_array

def test_swt_dtypes():
    wavelet = pywt.Wavelet('haar')
    for dt_in, dt_out in zip(dtypes_in, dtypes_out):
        errmsg = "wrong dtype returned for {0} input".format(dt_in)

        # swt
        x = np.ones(8, dtype=dt_in)
        (cA2, cD2), (cA1, cD1) = pywt.swt(x, wavelet, level=2)
        assert_(cA2.dtype == cD2.dtype == cA1.dtype == cD1.dtype == dt_out,
                "swt: " + errmsg)

        # swt2
        x = np.ones((8, 8), dtype=dt_in)
        cA, (cH, cV, cD) = pywt.swt2(x, wavelet, level=1)[0]
        assert_(cA.dtype == cH.dtype == cV.dtype == cD.dtype == dt_out,
                "swt2: " + errmsg)