def waverec_lin(x,wavelet,mode="per",sc_lengths=None):
    """Reconstructs the timeseries from the wavelet-coefficients as returned by wavedec_lin.""" 
    assert len(x.shape) in [1,2], "Only 1d and 2d-arrays supported up to now!"
    if sc_lengths != None:
        assert sum(sc_lengths) == x.shape[0], "The values given for sc_lengths don't make sense!"
    else:
        sc_lengths = []
        l = x.shape[0]
        while l>1:
            l/=2
            sc_lengths.insert(0,l)
        sc_lengths.insert(0,l)
        assert sum(sc_lengths) == x.shape[0], "The values automatically calculated for sc_lengths don't make sense!"
            
    if len(x.shape) == 1:
        wts = []
        for j in range(len(sc_lengths)):
            offset = sum(sc_lengths[:j])
            wts.append(x[offset:offset+sc_lengths[j]])
        rv = waverec(wts,wavelet,mode=mode)
    else: #len(x.shape)==2
        for i in range(x.shape[1]):
            wts = []
            for j in range(len(sc_lengths)):
                offset = sum(sc_lengths[:j])
                wts.append(x[offset:offset+sc_lengths[j],i])
            if i==0:
                ar = waverec(wts,wavelet,mode=mode)
                rv = n.zeros((ar.shape[0],x.shape[1]),"d")
                rv[:,i] = ar
            rv[:,i] = waverec(wts,wavelet,mode=mode)
    return rv

def wavelet_trasnform(sig):
    mode = pywt.Modes.smooth
    w = 'coif3'
    w = pywt.Wavelet(w)
    
#    hasil=[]
#    for n in range(data.shape[0]):
        #dd = data[n,1:24]
    ca = []
    cd = []
    for i in range(5):
            (a, d) = pywt.dwt(sig, w, mode)
            ca.append(a)
            cd.append(d)
    rec_a = []
    rec_d = []
    for i, coeff in enumerate(ca):
        coeff_list = [coeff, None] + [None] * i
        rec_a.append(pywt.waverec(coeff_list, w))
    for i, coeff in enumerate(cd):
        coeff_list = [None, coeff] + [None] * i
        rec_d.append(pywt.waverec(coeff_list, w))
#    hasil.append(rec_a[0])
        
    return rec_a[0];

def ApproximationsV(data, family, levels):
    """
    get approximation reconstrutions at different levels
    for a DWT family.
    returns levels+1 arrays with A[0]=full reconstruction,
    and A[1]=first approx, A[levels] is smoothest
    """
    # subtract off mean data
    meandata=np.mean(data)
    #meandata=0.0
    # get DWT coefficients
    coeffs = pywt.wavedec(data-meandata, family, mode='sym',level=Nlevels)
    lcoeffs=len(coeffs)
    for i,l in enumerate(coeffs):
        vl=np.var(l)
        l[:]=vl
        coeffs[i]=l
    #print "len coeffs",lcoeffs
    #for i in coeffs: print len(i)
    # reconstruct approximations
    A=[]
    c=pywt.waverec(coeffs,family,mode='sym')
    A.append(np.array(c)+meandata)
    for j in range(Nlevels,0,-1):
        coeffs[j][0:]=0.0
        c=pywt.waverec(coeffs,family,mode='sym')
        A.append(np.array(c)+meandata)
    return A

def wave_stein_filter(y, sigma, tau, w):
    coeffs = pywt.wavedec(y, w)
    threshold = sigma * tau
    hcoeffs = []
    for scale, x in enumerate(coeffs):
            hcoeffs.append(stein_thresholding(x, threshold))
    return pywt.waverec(hcoeffs, w)

def wave_semisoft_filter(y, sigma, tau, w, mu):
    coeffs = pywt.wavedec(y, w)
    threshold = sigma * tau
    hcoeffs = []
    for scale, x in enumerate(coeffs):
            hcoeffs.append(thresholding_semisoft(x, threshold, mu))
    return pywt.waverec(hcoeffs, w)

def haar(vals, p=80):
	hC = pywt.wavedec(vals,'haar')
	cutVal = np.percentile(np.abs(np.concatenate(hC)), p)
	for A in hC:
		A[np.abs(A) < cutVal] = 0
	tVals = pywt.waverec(hC,'haar')
	return tVals[:len(vals)]

def denoise(nblck,filename,mode='sym', wv='sym5' ):
    from statsmodels.robust import mad
    #noisy_coefs = pywt.wavedec(nblck, 'sym5', level=5, mode='per')
    noisy_coefs = pywt.wavedec(nblck, wavelet=wv,   mode=mode) #level=5,  #dwt is for single level decomposition; wavedecoding is for more levels
    sigma = mad(noisy_coefs[-1])
    #uthresh=np.std(ca)/2
    uthresh = sigma*np.sqrt(2*np.log(len(nblck)))
    denoised = noisy_coefs[:]
    denoised[1:] = [pywt.threshold(i, value=uthresh,mode='soft') for i in denoised[1:]]
    signal = pywt.waverec(denoised, wavelet=wv, mode=mode)
    from matplotlib import pyplot as plt
    fig, axes = plt.subplots(1, 2, sharey=True, sharex=True,figsize=(8,4))
    ax1, ax2 = axes
    
    ax1.plot(signal)
    #ax1.set_xlim(0,2**10)
    ax1.set_title("Recovered Signal")
    ax1.margins(.1)
    
    ax2.plot(nblck)
    ax2.set_title("Noisy Signal")
    
    for ax in fig.axes:
        ax.tick_params(labelbottom=False, top=False, bottom=False, left=False,       right=False)
    fig.tight_layout()
    fig.savefig(filename+'_'+wv+'.pdf')
    plt.clf()
    return signal

def generateCompression(cumHist, wavelet, compression):
	# Do a full multilevel wavelet decomposition until some level/resolution
	coeffs = pywt.wavedec(cumHist, wavelet)

	# Reduce trailing equal values from vectors
	vectors, lengths = reduceVectors(coeffs, 0.01)

	# Overwrite least significant detail coefficients in coeffs with 0-arrays (simulating removal)
	if compression < len(vectors):
		for i in range(compression,len(vectors)):
			arraySize = len(vectors[i])
			vectors[i] = np.zeros(arraySize)
	
	# Generate statistics
	stored = len(lengths)
	traildel = 0
	for i in range(0,min(len(vectors), compression)):
		stored += len(vectors[i])
		traildel += lengths[i] - len(vectors[i])
	
	# Restore trailing equal values
	restoredCoeffs = restoreVectors(vectors, lengths)

	# Do a full multilevel wavelet reconstruction from coeffs (zero-padded)
	recCumHist = pywt.waverec(restoredCoeffs, wavelet)
	return restoredCoeffs, recCumHist, stored, traildel, coeffs

	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 _call(self, coeff):
        """Compute the discrete 1D, 2D or 3D inverse wavelet transform.

        Parameters
        ----------
        coeff : `DiscreteLpVector`

        Returns
        -------
        arr : `DiscreteLpVector`

        """
        if len(self.range.shape) == 1:
            coeff_list = array_to_pywt_coeff(coeff, self.size_list)
            x = pywt.waverec(coeff_list, self.wbasis, self.mode)
            return self.range.element(x)
        elif len(self.range.shape) == 2:
            coeff_list = array_to_pywt_coeff(coeff, self.size_list)
            x = pywt.waverec2(coeff_list, self.wbasis, self.mode)
            return self.range.element(x)
        elif len(self.range.shape) == 3:
            coeff_dict = array_to_pywt_coeff(coeff, self.size_list)
            x = wavelet_reconstruction3d(coeff_dict, self.wbasis, self.mode,
                                         self.nscales)
            return self.range.element(x)

def soft_threshold(bins, tau, level=1, plot=False):
        coeffs = pywt.wavedec(bins, 'haar', level=level)
        old_coeffs = copy.deepcopy(coeffs) if plot else None
        coeffs[1:] = [sign(a) * np.maximum(0, abs(a)-tau) for a in coeffs[1:] ]
        filtered = pywt.waverec(coeffs, 'haar')

        if plot:
                from matplotlib import pyplot as pl
                fig = pl.figure(figsize=(15,8))
                fig.suptitle("Level %d Wavelet Decomposition" % level)
                fig.subplots_adjust(hspace=0.3, wspace=0.3)

                ax = fig.add_subplot(2,1,1)
                ax.plot(bins, label='Raw')
                ax.plot(filtered, label='Filtered')
                ax.legend()

                ax = fig.add_subplot(2,level,1+level)
                ax.set_title("Approximation coefficients")
                ax.plot(coeffs[0])
                for l in range(1,level):
                        ax = fig.add_subplot(2, level, 1+l+level)
                        ax.plot(old_coeffs[l])
                        ax.axhline(-tau, color='g'); ax.axhline(+tau, color='g')
                        ax.set_title("Level %d detail coefficients" % l)

                fig.savefig('wavelet-analysis.pdf')

        return filtered

    def filter(self):

        if self.level > self.max_dec_level():
            clevel = self.max_dec_level()
        else:
            clevel = self.level

        # decompose
        coeffs = pywt.wavedec(self.sig, pywt.Wavelet(self.wt), \
                              mode=self.mode, \
                              level=clevel)

        # threshold evaluation
        th = sqrt(2 * log(len(self.sig)) * power(self.sigma, 2))

        # thresholding
        for (i, cAD) in enumerate(coeffs):
            if i == 0:
                continue
            coeffs[i] = sign(cAD) * pywt.thresholding.less(abs(cAD), th)

        # reconstruct
        rec_sig = pywt.waverec(coeffs, pywt.Wavelet(self.wt), mode=self.mode)
        if len(rec_sig) == (len(self.sig) + 1):
            self.sig = rec_sig[:-1]

 def denoise(self, data, wavelet):
     noiseSigma = median(absolute(data - median(data))) / 0.6745
     levels = int(floor(log(len(data))))
     WC = pywt.wavedec(data, wavelet, level=levels)
     threshold = noiseSigma * sqrt(2 * log(len(data)))
     NWC = map(lambda x: pywt.thresholding.hard(x, threshold), WC)
     return pywt.waverec(NWC, wavelet)

    def _wm_reverse(self, wd):
        if __debug__:
            debug('MAP', "Performing iDWT")
        signal = None
        wd_offsets = [0] + list(np.cumsum(self.lengths))
        nlevels = len(self.lengths)
        Ntime_points = self._intimepoints #len(time_points)
        # unfortunately sometimes due to padding iDWT would return longer
        # sequences, thus we just limit to the right ones

        for indexes in _get_indexes(wd.shape, self._dim):
            if __debug__:
                debug('MAP_', " %s" % (indexes,), lf=False, cr=True)
            wd_sample = wd[indexes]
            wd_coeffs = [wd_sample[wd_offsets[i]:wd_offsets[i+1]] for i in xrange(nlevels)]
            # need to compose original list
            time_points = pywt.waverec(
                wd_coeffs, wavelet=self._wavelet, mode=self._mode)
            if signal is None:
                newdim = list(wd.shape)
                newdim[self._dim] = Ntime_points
                signal = np.empty(newdim)
            signal[indexes] = time_points[:Ntime_points]
        if __debug__:
            debug('MAP_', "")
            debug('MAP', "Done iDWT. Total size %s" % (signal.shape, ))
        return signal

def test_waverec_axis_db2():
    # test for fix to issue gh-293
    rstate = np.random.RandomState(0)
    data = rstate.standard_normal((16, 16))
    for axis in [0, 1]:
        coefs = pywt.wavedec(data, 'db2', axis=axis)
        rec = pywt.waverec(coefs, 'db2', axis=axis)
        assert_allclose(rec, data, atol=1e-14)

def visu_shrink_denoise(f, wave, v):
    true_coefs = pywt.wavedec(f, wave, level=2)
    res = [visu_shrink(coef) for coef in true_coefs]
    den_coefs, thresholds = zip(*res)
    den_coefs = list(den_coefs)
    #den_coefs.insert(0, true_coefs[0])
    f_denoised = pywt.waverec(den_coefs, wave)
    return f_denoised

 def rmatvec(x):
     x_list = []
     count = 0
     for el in b:
         n_el = np.asarray(el).size
         x_list.append(np.array(x[count:count+n_el]))
         count += n_el
     return pywt.waverec(x_list, wavelet, mode=mode)

def sure_shrink_denoise(f, wave, v):
    true_coefs = pywt.wavedec(f, wave, level=None)
    res = [sure_shrink_fast(coef) for coef in true_coefs[1:]]
    den_coefs, thresholds = zip(*res)
    den_coefs = list(den_coefs)
    den_coefs.insert(0, true_coefs[0])
    f_denoised = pywt.waverec(den_coefs, wave)
    return f_denoised

def idwt(E):
	# print E
	E_rep=E.tolist()[0:2]
	for i in range(1,int(pow(len(E),0.5))):
		E_rep.append(E[pow(2,i):pow(2,i+1)])
	E_idwt=pywt.waverec(E_rep,'db1')
	# print E_idwt
	return matrix(E_idwt)

def waveletDenoise(data, threshold=1, wavelet='haar'):
    """ Wavelet denoise test """
    # Decompose
    wl_dec = pywt.wavedec(data, pywt.Wavelet(wavelet))
    # Threshold
    wl_thr = map(lambda x: pywt.thresholding.soft(x, threshold), wl_dec)
    # Reconstruct
    return pywt.waverec(wl_thr, wavelet)