def wavedec_lin(x,wavelet,mode="per",level=None,return_sl=False): 
    """Performs a wavelet-decomposition using pywt.wavedec, but returns the coefficients
    as one 1d-array.
    Uses different default-values for mode and level.
    Behaviour for the case "len(x) not a power of 2" is not complete!"""       
    assert len(x.shape) in [1,2], "Only 1d and 2d-arrays supported up to now!"
    if level == None:
        level = int(round(n.log2(x.shape[0])))
    if len(x.shape) == 1:
        wts = wavedec(x,wavelet,mode=mode,level=level)
        sc_lengths = [len(x) for x in wts]
        rv = n.zeros((sum(sc_lengths)),"d")
        for j in range(len(wts)):
            offset = sum(sc_lengths[:j])
            rv[offset:offset+sc_lengths[j]]=wts[j]
    else: #len(x.shape)==2
        wts = wavedec(x[:,0],wavelet,mode=mode,level=level)
        sc_lengths = [len(i) for i in wts]
        rv = n.zeros((sum(sc_lengths),x.shape[1]),"d")
        for i in range(x.shape[1]):
            if i>0:
                wts = wavedec(x[:,i],wavelet,mode=mode,level=level)
            for j in range(len(wts)):
                offset = sum(sc_lengths[:j])
                rv[offset:offset+sc_lengths[j],i]=wts[j]
    if return_sl:
        return rv, sc_lengths
    else:
        return rv

    def fit(self):
        '''
        First finds the wavelet in types that fits the data the best => The smallest Euclidean distance between the 
        approximation signal and under-sampled signal
        After finding the best Wavelet, uses it to decompose data into levels detail signals and builds a AR(p)
        model per approximation signal and detail signals using order as p for each corresponding  
        '''
        self.bestType = 'db1'
        bestDist = np.inf

        for t in self.types:
            approx_levels = pywt.wavedec(self.data, wavelet=t, level=self.levels)
            idx = np.int_(np.linspace(0, self.data.shape[0]-1, num=len(approx_levels[0])))
            samples = self.data[idx]
            dist = np.linalg.norm(approx_levels[0]-samples)
            if dist < bestDist:
                bestDist = dist
                self.bestType = t
                
        self.coefs = pywt.wavedec(self.data, wavelet=self.bestType, level=self.levels)
        for i in range(len(self.order)):
#             model = AR_model(approx_levels[i], order=self.order[i])
            model = Markov_model(approx_levels[i])
            model.fit()
            self.models.append(model)

def _weighted_retrieve(data, genome, loci, prediction_steps, spinup, weight_func):
    l = loci
    (temps, loads) = (data['Temperature'], data['Load'])

    idx = _grow_tree(data, genome, loci, prediction_steps)
    window = genome[l.hindsight]
    test_starts, test_ends = _get_test_period(data)
    
    query_loads_norm = sg.utils.Normalizer(loads[test_starts-window:test_starts])
    query_loads = [item for sublist in pywt.wavedec(query_loads_norm.normalized, 'haar')
                   for item in sublist ][:idx.properties.dimension/2]

    query_temps_norm = sg.utils.Normalizer(temps[test_starts-window:test_starts])
    query_temps = [item for sublist in pywt.wavedec(query_temps_norm.normalized, 'haar')
                   for item in sublist ][:idx.properties.dimension/2]

    query = weight_func(query_loads) + weight_func(query_temps)
    # If we find matches where there is a gap in the timeseries (which will throw an exception), we look for the next best match.
    num_matches = 1
    while True:
        match_date = list(idx.nearest(tuple(query), num_matches, objects="raw"))[-1]
        end_date = match_date + dt(hours=genome[l.hindsight]+prediction_steps - 1)
        period = sg.utils.Normalizer(loads[match_date:end_date]).normalized
        prediction = query_loads_norm.expand(period[-prediction_steps:])
        try:
            result = pd.TimeSeries(data=prediction.values,
                                   index=data[test_starts:test_ends+1].index)
            idx.close()
            return result
        except:
            num_matches += 1
            print 'Time gap encountered, we will try match number', num_matches

 def save_samples(self):
     """
     Save all the generated samples in a samples.txt file
     """
     
     doc_count = sim_count = 1
     with open("%ssets/learn_test_set.txt" % settings.MEDIA_ROOT, "w") as lt_file:
         for sample in self.learn_set:
             if doc_count % 9 == 0:
                 with open("%ssets/sim_sample%s.txt" % (settings.MEDIA_ROOT, sim_count), "w") as sim_doc:
                     self._write_sample(sim_doc, sample)
                 sim_count += 1
                 self.targets.pop(doc_count-1)
             else:
                 sample = wavedec(sample, self.wavelet, level=5)
                 sample = [i/10 for i in list(sample[0])]
                 self._write_sample(lt_file, sample)
                 lt_file.write(";")
             doc_count += 1
         for sample in self.test_set:
             sample = wavedec(sample, self.wavelet, level=5)
             sample = [i/10 for i in list(sample[0])]
             self._write_sample(lt_file, sample)
             lt_file.write(";")
     with open("%ssets/target_set.txt" % settings.MEDIA_ROOT, "w") as target_doc:
         self._write_sample(target_doc, self.targets)

def decomposeField(dataframe, fieldName, groupFieldName, maxCoef) :
    coeffs = {} #Coefficients for each group
    maxLen = 0  #Largest list of coefficients

    #Grouped Case
    try:
        grouped = dataframe.groupby(groupFieldName)
        for name, group in grouped:
            #Collect coefficients for each group
            coeffs[name] = pywt.wavedec(group[fieldName], 'db1', level=2)[0].tolist()[:maxCoef]
            maxLen = max(maxLen,len(coeffs[name]))
            #Non-grouped case
    except KeyError:
        #No group.  One row of coefficients
        coeffs[0] = pywt.wavedec(dataframe[fieldName], 'db1', level=2)[0].tolist()[:maxCoef]
        maxLen = len(coeffs[0])

    # Ensures all rows of coefficients are the same length.  
    # Populates anything shorter with nan
    for coef in coeffs:
        coeffs[coef] = coeffs[coef] + [float('nan')]*(maxLen-len(coeffs[coef]))
    #Assign names of coefficients using the original field name as a prefix
    names = [fieldName + str(i) for i in range(maxLen)]    #note change from original

    #Transpose & return
    coeffD = pd.DataFrame(coeffs)
    coeffT = coeffD.T
    coeffT.columns = names
    return coeffT

def wavepower_db_secondhalf(x,wavelet,mode="per",level=None):
    """Unused? Maybe remove it."""
    if level == None:
        level = int(round(n.log2(x.shape[0])))
    if len(x.shape) == 1:
        wts = wavedec(x,wavelet,mode=mode,level=level)
        sc_lengths = [len(x) for x in wts]
        rv = n.zeros((sum(sc_lengths)/2),"d")
        for j in range(len(wts)):
            wts_power = wts[j]**2
            idx_norm = max( 0 , int(wts_power.shape[0]/2)-1 )
                #print "idx_norm:", idx_norm
            wts_power /= wts_power[idx_norm]
            wts_power = 10* np.log(wts_power)
            offset = sum(sc_lengths[:j])/2
            print rv[offset:offset+sc_lengths[j]/2].shape, wts_power[idx_norm+1:].shape
            rv[offset:offset+sc_lengths[j]/2]=wts_power[idx_norm+1:]
    else: #len(x.shape)==2
        return NotImplemented
        wts = wavedec(x[:,0],wavelet,mode=mode,level=level)
        sc_lengths = [len(i) for i in wts]
        rv = n.zeros((sum(sc_lengths),x.shape[1]),"d")
        for i in range(x.shape[1]):
            if i>0:
                wts = wavedec(x[:,i],wavelet,mode=mode,level=level)
            for j in range(len(wts)):
                wts_power = wts[j]**2
                if normalise:
                    idx_norm = max( 0 , int(wts_power.shape[0]*norm_fraction)-1 )
                    wts_power /= wts_power[idx_norm]
                offset = sum(sc_lengths[:j])
                rv[offset:offset+sc_lengths[j],i]=wts_power[:]
    return rv

 def wave_ica(x):
     ica = decomposition.FastICA(max_iter=10000)
     # calculate wavelet transform of each
     w1 = pywt.wavedec(x[:, 0], 'haar', level=12)
     w2 = pywt.wavedec(x[:, 1], 'haar', level=12)
     # calculate ica between components
     t = [ica.fit_transform(np.array([lev1, lev2]).T)
          for lev1, lev2 in zip(w1, w2)]
     return np.array(t)

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 extract_all_wcs_by_maxchan(self, wavelet='haar'):
        """Extract wavelet coefficients from all spikes, store them as spike attribs.
        Find optimum coeffs for each chan, then average across all chans to find
        globally optimum coeffs"""
        # TODO: add multiprocessing
        nkeep = 5 # num of top wavelet coeffs to keep
        sort = self.sort
        spikes = sort.spikes # struct array
        wavedata = sort.wavedata
        nspikes = len(spikes)
        #ncoeffs = 53 # TODO: this only applies for V of length 50, stop hardcoding
        #ncoeffs = len(self.ksis)
        nt = wavedata.shape[2]
        ncoeffs = len(np.concatenate(pywt.wavedec(wavedata[0, 0], wavelet)))

        wcs = {}
        maxchans = np.unique(spikes['chan'])
        nmaxchans = len(maxchans)
        for maxchan in maxchans:
            wcs[maxchan] = [] # init dict of lists, indexed by spike maxchan
        flatwcs = np.zeros((nspikes, ncoeffs))

        t0 = time.time()
        for spike, wd in zip(spikes, wavedata):
            nchans = spike['nchans']
            chans = spike['chans'][:nchans]
            maxchan = spike['chan']
            maxchani = int(np.where(chans == maxchan)[0])
            #chanis = det.chans.searchsorted(chans) # det.chans are always sorted
            #wd = wd[:nchans] # unnecessary?
            V = wd[maxchani]
            coeffs = np.concatenate(pywt.wavedec(V, wavelet)) # flat array of wavelet coeffs
            wcs[maxchan].append(coeffs)
            flatwcs[spike['id']] = coeffs
        ks = np.zeros((nmaxchans, ncoeffs))
        p = np.zeros((nmaxchans, ncoeffs))
        for maxchani, maxchan in enumerate(maxchans):
            wcs[maxchan] = np.asarray(wcs[maxchan])
            for i in range(ncoeffs):
                ks[maxchani, i], p[maxchani, i] = scipy.stats.kstest(wcs[maxchan][:, i], 'norm')
        ## TODO: weight the KS value from each maxchan according to the nspikes for that
        ## maxchan!!!!!
        ks = ks.mean(axis=0)
        p = p.mean(axis=0)
        ksis = ks.argsort()[::-1] # ks indices sorted from biggest to smallest ks values
        # assign as params in spikes struct array
        for coeffi in range(nkeep): # assign first nkeep
            spikes['w%d' % coeffi] = flatwcs[:, ksis[coeffi]]
        print("Extracting wavelet coefficients from all %d spikes took %.3f sec" %
             (nspikes, time.time()-t0))
        return wcs, flatwcs, ks, ksis, p

        def dot(self, mat):
            m = []

            if mat.shape[0] != mat.size:
                for i in xrange(mat.shape[1]):
                    c = pywt.wavedec(mat[:, i], self.name, level=self.level)
                    self.sizes.append(map(len, c))
                    c = np.concatenate(c)
                    m.append(c)
                return np.asarray(m).T
            else:
                c = pywt.wavedec(mat, self.name, level=self.level)
                self.sizes.append(map(len, c))
                return np.concatenate(c)

def dwt_eeg_video(video_eeg_data, electrode_count, electrode_indexes):
    """
    Use Discrete wavelet transform (DWT) to compute alpha and beta band of signal. Compute power of alpha and beta band
    and also valence and arousal values.
    :param video_eeg_data: eeg signal
    :param electrode_count: number of electrodes
    :param electrode_indexes: indexes of electrodes, usually just range(0, electrode_count)
    :return: array of floats, shape [electrode_count*2 + 2, 1]
             power of alpha and beta band of individual electrodes, valence and arousal values computed from eeg signal
    notes: this function should be split into more in the future
    """

    data_final = np.empty(electrode_count*2 + 2)

    alphaArray = []
    betaArray = []
    counter = 0
    for electrodeIndex in electrode_indexes:
        coeffs = pywt.wavedec(video_eeg_data[electrodeIndex], 'db2', level=3)
        a3, d3, d2, d1 = coeffs

        coeffs = pywt.wavedec(d3, 'db2', level=1)

        alpha, beta = coeffs
        alphaArray.append(power_of_signal(alpha))
        data_final[counter] = power_of_signal(alpha)

        beta = pywt.idwt(d2,sig.resample(beta,d2.__len__()),'db2')
        betaArray.append(power_of_signal(beta))
        data_final[counter+1] = power_of_signal(beta)

        counter += 2

    F3alpha = alphaArray[0]
    F4alpha = alphaArray[1]
    AF3alpha = alphaArray[2]
    AF4alpha = alphaArray[3]
    F3beta = betaArray[0]
    F4beta = betaArray[1]
    AF3beta = betaArray[2]
    AF4beta = betaArray[3]

    valence = (F4alpha/F4beta) - (F3alpha/F3beta)
    arousal = (F3beta+F4beta+AF3beta+AF4beta) / (F3alpha+F4alpha+AF3alpha+AF4alpha)

    data_final[counter] = valence
    data_final[counter+1] = arousal

    return data_final

    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 WT(data, wavelet, mode='sym'):
    """Perfroms a batch 1D wavelet transform.

    Parameters
    ----------
    data : array
        (n_vars, n_obs, n_contacts) array where `n_vars` is the number of
        variables (vector dimensions), `n_obs` the number of observations
        and `n_contacts` is the number of contacts. Only 3D arrays are
        accepted.
    wavelet : string or pywt.Wavelet
        wavelet to be used to perform the transform
    mode : string, optional
        signal extension mode (see modes in PyWavelets documentation)

    Returns
    -------
    data : array
        (n_coeffs, n_obs, n_contacts) 1D wavelet transform of each vector
        of the input data array. `pywt.wavedec` is used to perform the
        transform. For every vector of the input array, a 1D transformation
        is returned of the form [cA_n, cD_n, cD_n-1, ..., cD_n2, cD_n1]
        where cA_n and cD_n are approximation and detailed coefficients of
        level n. cA_n and cD_n's are stacked together in a single vector.

    Notes
    -----

    PyWavelets documentation contains more detailed information on the
    wavelet transform.

    """
    # TODO: complete docstring, add wavelet type checking,
    # think about dependencies

    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

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

    n_samples = data.shape[0]
    n_spikes = data.shape[1]
    n_contacts = data.shape[2]
    n_features = full_coeff_len(n_samples, wavelet.dec_len, mode)
    new_data = np.empty((n_features, n_spikes, n_contacts))

    for i in xrange(n_spikes):
        for c in xrange(n_contacts):
            coeffs = wt.wavedec(data[:, i, c], wavelet, mode)
            new_data[:, i, c] = np.hstack(coeffs)

    return new_data

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_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 dwtTransform(x):
    xdwt = []
    coeffs = pywt.wavedec(x,'haar')
    for coeff in  coeffs:
        for a in coeff:
            xdwt.append(a)
    return xdwt

    def calculate_mra(self, wavelet='db10', mode='per'):
        """
        Creates an MRA wavelet tree on the recording.
        
        Args:
            wavelet (str): wavelet to use. Any string supported by PyWavelets will work.
            mode (str): method for handling overrun. Default "per," start over at the beginning of the waveform
                (periodic).
        """
        
        self.wavelet, self.mode = wavelet, mode
        self.dwt = pywt.wavedec(self.wav, wavelet, mode=mode, level=int(np.log2(len(self.wav))) + 1)
        
        self.root = None
        self.nodes = []
        self.wavelet = wavelet
        self.mode = mode
        parents = [None]

        for i in range(len(self.dwt)):
            nodes = []

            for j in range(len(self.dwt[i])):
                if j > 0:
                    nodes.append(MRANode(self.dwt[i][j], parents[j / 2], nodes[j - 1], j % 2, i))
                else:
                    nodes.append(MRANode(self.dwt[i][j], parents[j / 2], None, j % 2, i))

            nodes[0].predecessor = nodes[-1]
            parents = nodes
            self.nodes.extend(nodes)
            if i is 0: self.root = nodes[0]

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 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