<commit_before># -*- coding: utf-8 -*-
"""
Created on Mon Oct  7 13:22:58 2013

@author: hari 
"""
import nitime.algorithms as alg
import numpy as np
from math import ceil
import scipy as sci
from scipy import linalg
from .utils import logger, verbose

@verbose
def mtplv(x, params, verbose = None):
    """Multitaper Phase-Locking Value
    
    Parameters
    ----------
    x - NumPy Array 
        Input Data (channel x trial x time) or (trials x time)
    params - Dictionary of parameter settings
      params['Fs'] - sampling rate
      
      params['tapers'] - [TW, Number of tapers]
      
      params['fpass'] - Freqency range of interest, e.g. [5, 1000]
      
      params['pad'] - 1 or 0, to pad to the next power of 2 or not
      
      params['itc'] - 1 for ITC, 0 for PLV
    verbose : bool, str, int, or None
        The verbosity of messages to print. If a str, it can be either DEBUG,
        INFO, WARNING, ERROR, or CRITICAL.
        
    Returns
    -------
    (plvtap, f): Tuple
        plvtap - Multitapered phase-locking estimate (channel x frequency)
        
        f - Frequency vector matching plvtap
    """
    logger.info('Running Multitaper PLV Estimation')
    if(len(x.shape) == 3):
        timedim = 2
        trialdim = 1
        nchans = x.shape[0]
        ntrials = x.shape[trialdim]
        logger.info('The data is of format %d channels x %d trials x time',
                    nchans, ntrials)
    elif(len(x.shape) == 2):
        timedim = 1
        trialdim = 0
        ntrials = x.shape[trialdim]
        nchans = 1
        logger.info('The data is of format %d trials x time (single channel)',
                    ntrials)
    else:
        logger.error('Sorry, The data should be a 2 or 3 dimensional array')
        
    # Calculate the tapers
    ntaps = params['tapers'][1]
    TW = params['tapers'][0]
    w,conc = alg.dpss_windows(x.shape[timedim],TW,ntaps)
    
    # Make space for the PLV result
    Fs = params['Fs']
    nfft = int(2**ceil(sci.log2(x.shape[timedim])))
    f = np.arange(0.0,nfft,1.0)*Fs/nfft
    plvtap = np.zeros((ntaps,nchans,nfft))
    
    
    for k,tap in enumerate(w):
        logger.info('Doing Taper #%d', k)
        xw = sci.fft(tap*x,n = nfft, axis = timedim)
        
        if(params['itc'] == 0):
            plvtap[k,:,:] = abs((xw/abs(xw)).mean(axis = trialdim))**2
        else:
            plvtap[k,:,:] = ((abs(xw.mean(axis = trialdim))**2)/
            ((abs(xw)**2).mean(axis = trialdim)))
            
    plvtap = plvtap.mean(axis = 0).squeeze()
    ind = (f > params['fpass'][0]) & (f < params['fpass'][1])
    plvtap = plvtap[:,ind]
    f = f[ind]
    return (plvtap,f)

@verbose
def mtspec(x,params, verbose = None):
    """Multitaper Spectrum and SNR estimate
    
    Parameters
    ----------
    x - NumPy Array
        Input data (channel x trial x time) or (trials x time)
    
    params - Dictionary of parameter settings
      params['Fs'] - sampling rate
      
      params['tapers'] - [TW, Number of tapers]
      
      params['fpass'] - Freqency range of interest, e.g. [5, 1000]
      
      params['pad'] - 1 or 0, to pad to the next power of 2 or not
    verbose : bool, str, int, or None
        The verbosity of messages to print. If a str, it can be either DEBUG,
        INFO, WARNING, ERROR, or CRITICAL.
      
    Returns
    -------
    (S, N ,f): Tuple
        S - Multitapered spectrum (channel x frequency)
        
        N - Noise floor estimate
        
        f - Frequency vector matching S and N
    """
    
    logger.info('Running Multitaper Spectrum and Noise-floor Estimation')
    if(len(x.shape) == 3):
        timedim = 2
        trialdim = 1
        ntrials = x.shape[trialdim]
        nchans = x.shape[0]
        logger.info('The data is of format #d channels x %d trials x time',
                     nchans, ntrials)
    elif(len(x.shape) == 2):
        timedim = 1
        trialdim = 0
        ntrials = x.shape[trialdim]
        nchans = 1
        logger.info('The data is of format %d trials x time (single channel)',
                    ntrials)
    else:
        logger.error('Sorry! The data should be a 2 or 3 dimensional array!')
        
    # Calculate the tapers
    ntaps = params['tapers'][1]
    TW = params['tapers'][0]
    w,conc = alg.dpss_windows(x.shape[timedim],TW,ntaps)
    
    # Make space for the results
    Fs = params['Fs']
    nfft = int(2**ceil(sci.log2(x.shape[timedim])))
    f = np.arange(0.0,nfft,1.0)*Fs/nfft
    S = np.zeros((ntaps,nchans,nfft))
    N = np.zeros((ntaps,nchans,nfft))
    
    
    for k,tap in enumerate(w):
        logger.info('Doing Taper #%d',k)
        xw = sci.fft(tap*x,n = nfft, axis = timedim)
        # randph = sci.rand(nchans,ntrials,nfft)*2*sci.pi
        randsign = np.ones((nchans,ntrials,nfft))
        randsign[:,np.arange(0,ntrials,2),:] = -1
        S[k,:,:] = abs(xw.mean(axis = trialdim))
        # N[k,:,:] = abs((xw*sci.exp(1j*randph)).mean(axis = trialdim))
        N[k,:,:] = abs((xw*(randsign.squeeze())).mean(axis = trialdim))
            
    # Average over tapers and squeeze to pretty shapes        
    S = S.mean(axis = 0).squeeze() 
    N = N.mean(axis = 0).squeeze()
    ind = (f > params['fpass'][0]) & (f < params['fpass'][1])
    S = S[:,ind]
    N = N[:,ind]
    f = f[ind]
    return (S,N,f)

@verbose      
def mtcpca(x,params, verbose = None):
    """Multitaper complex PCA and PLV
    
    Parameters
    ----------
    x - NumPy Array
        Input data (channel x trial x time)
    
    params - Dictionary of parameter settings
      params['Fs'] - sampling rate
      <commit_msg>ENH: Allowing different types of noise-floor estimates<commit_after># -*- coding: utf-8 -*-
"""
Created on Mon Oct  7 13:22:58 2013

@author: hari 
"""
import nitime.algorithms as alg
import numpy as np
from math import ceil
import scipy as sci
from scipy import linalg
from .utils import logger, verbose

@verbose
def mtplv(x, params, verbose = None):
    """Multitaper Phase-Locking Value
    
    Parameters
    ----------
    x - NumPy Array 
        Input Data (channel x trial x time) or (trials x time)
    params - Dictionary of parameter settings
      params['Fs'] - sampling rate
      
      params['tapers'] - [TW, Number of tapers]
      
      params['fpass'] - Freqency range of interest, e.g. [5, 1000]
      
      params['pad'] - 1 or 0, to pad to the next power of 2 or not
      
      params['itc'] - 1 for ITC, 0 for PLV
    verbose : bool, str, int, or None
        The verbosity of messages to print. If a str, it can be either DEBUG,
        INFO, WARNING, ERROR, or CRITICAL.
        
    Returns
    -------
    (plvtap, f): Tuple
        plvtap - Multitapered phase-locking estimate (channel x frequency)
        
        f - Frequency vector matching plvtap
    """
    logger.info('Running Multitaper PLV Estimation')
    if(len(x.shape) == 3):
        timedim = 2
        trialdim = 1
        nchans = x.shape[0]
        ntrials = x.shape[trialdim]
        logger.info('The data is of format %d channels x %d trials x time',
                    nchans, ntrials)
    elif(len(x.shape) == 2):
        timedim = 1
        trialdim = 0
        ntrials = x.shape[trialdim]
        nchans = 1
        logger.info('The data is of format %d trials x time (single channel)',
                    ntrials)
    else:
        logger.error('Sorry, The data should be a 2 or 3 dimensional array')
        
    # Calculate the tapers
    ntaps = params['tapers'][1]
    TW = params['tapers'][0]
    w,conc = alg.dpss_windows(x.shape[timedim],TW,ntaps)
    
    # Make space for the PLV result
    Fs = params['Fs']
    nfft = int(2**ceil(sci.log2(x.shape[timedim])))
    f = np.arange(0.0,nfft,1.0)*Fs/nfft
    plvtap = np.zeros((ntaps,nchans,nfft))
    
    
    for k,tap in enumerate(w):
        logger.info('Doing Taper #%d', k)
        xw = sci.fft(tap*x,n = nfft, axis = timedim)
        
        if(params['itc'] == 0):
            plvtap[k,:,:] = abs((xw/abs(xw)).mean(axis = trialdim))**2
        else:
            plvtap[k,:,:] = ((abs(xw.mean(axis = trialdim))**2)/
            ((abs(xw)**2).mean(axis = trialdim)))
            
    plvtap = plvtap.mean(axis = 0).squeeze()
    ind = (f > params['fpass'][0]) & (f < params['fpass'][1])
    plvtap = plvtap[:,ind]
    f = f[ind]
    return (plvtap,f)

@verbose
def mtspec(x,params, verbose = None):
    """Multitaper Spectrum and SNR estimate
    
    Parameters
    ----------
    x - NumPy Array
        Input data (channel x trial x time) or (trials x time)
    
    params - Dictionary of parameter settings
      params['Fs'] - sampling rate
      
      params['tapers'] - [TW, Number of tapers]
      
      params['fpass'] - Freqency range of interest, e.g. [5, 1000]
      
      params['pad'] - 1 or 0, to pad to the next power of 2 or not
      
      params['noisefloortype'] - (optional) 1: random phase, 
      0 (default): flip-phase
      
    verbose : bool, str, int, or None
        The verbosity of messages to print. If a str, it can be either DEBUG,
        INFO, WARNING, ERROR, or CRITICAL.
      
    Returns
    -------
    (S, N ,f): Tuple
        S - Multitapered spectrum (channel x frequency)
        
        N - Noise floor estimate
        
        f - Frequency vector matching S and N
    """
    
    logger.info('Running Multitaper Spectrum and Noise-floor Estimation')
    if(len(x.shape) == 3):
        timedim = 2
        trialdim = 1
        ntrials = x.shape[trialdim]
        nchans = x.shape[0]
        logger.info('The data is of format %d channels x %d trials x time',
                     nchans, ntrials)
    elif(len(x.shape) == 2):
        timedim = 1
        trialdim = 0
        ntrials = x.shape[trialdim]
        nchans = 1
        logger.info('The data is of format %d trials x time (single channel)',
                    ntrials)
    else:
        logger.error('Sorry! The data should be a 2 or 3 dimensional array!')
        
    # Calculate the tapers
    ntaps = params['tapers'][1]
    TW = params['tapers'][0]
    w,conc = alg.dpss_windows(x.shape[timedim],TW,ntaps)
    
    # Make space for the results
    Fs = params['Fs']
    nfft = int(2**ceil(sci.log2(x.shape[timedim])))
    f = np.arange(0.0,nfft,1.0)*Fs/nfft
    S = np.zeros((ntaps,nchans,nfft))
    N = np.zeros((ntaps,nchans,nfft))
    
    
    for k,tap in enumerate(w):
        logger.info('Doing Taper #%d',k)
        xw = sci.fft(tap*x,n = nfft, axis = timedim)
        
        S[k,:,:] = abs(xw.mean(axis = trialdim))
            
        if ('noisefloortype' in params) and (params['noisefloortype'] == 1):
            randph = sci.rand(nchans,ntrials,nfft)*2*sci.pi
            N[k,:,:] = abs((xw*sci.exp(1j*randph)).mean(axis = trialdim))
        else:
            randsign = np.ones((nchans,ntrials,nfft))
            randsign[:,np.arange(0,ntrials,2),:] = -1
            N[k,:,:] = abs((xw*(randsign.squeeze())).mean(axis = trialdim))
            
    # Average over tapers and squeeze to pretty shapes        
    S = S.mean(axis = 0).squeeze() 
    N = N.mean(axis = 0).squeeze()
    ind = (f > params['fpass'][0]) & (f < params['fpass'][1])
    S = S[:,ind]
    N = N[:,ind]
    f = f[ind]
    return (S,N,f)

@verbose      
def mtcpca(x,params, verbose = None):
    """Multitaper complex PCA and PLV
    
    Parameters
    ----------
    x - NumPy Array
        Input data (channel x trial x time)
    
    params - Dictionary of parameter settings
      params['Fs'] - sampling rate
      
{'fraction_non_alphanumeric': 0.13702548101924078, 'fraction_numerical': 0.014473912289824926, 'mean_word_length': 2.019890081130594, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '65896', 'n_tokens_mistral': 4198, 'n_tokens_neox': 3972, 'n_words': 1112}