#NEWLINE# def adaptive_weights( yk , eigvals , sides = 'onesided' , max_iter = 150 ) : #NEWLINE# #INDENT# '\n    Perform an iterative procedure to find the optimal weights for K\n    direct spectral estimators of DPSS tapered signals.\n\n    Parameters\n    ----------\n\n    yk : ndarray (K, N)\n       The K DFTs of the tapered sequences\n    eigvals : ndarray, length-K\n       The eigenvalues of the DPSS tapers\n    sides : str\n       Whether to compute weights on a one-sided or two-sided spectrum\n    max_iter : int\n       Maximum number of iterations for weight computation\n\n    Returns\n    -------\n\n    weights, nu\n\n       The weights (array like sdfs), and the\n       "equivalent degrees of freedom" (array length-L)\n\n    Notes\n    -----\n\n    The weights to use for making the multitaper estimate, such that\n    :math:`S_{mt} = \\sum_{k} |w_k|^2S_k^{mt} / \\sum_{k} |w_k|^2`\n\n    If there are less than 3 tapers, then the adaptive weights are not\n    found. The square root of the eigenvalues are returned as weights,\n    and the degrees of freedom are 2*K\n\n    ' #NEWLINE# from nitime.algorithms import mtm_cross_spectrum #NEWLINE# K = len ( eigvals ) #NEWLINE# if ( len ( eigvals ) < 3 ) : #NEWLINE# #INDENT# print ( '\n        Warning--not adaptively combining the spectral estimators\n        due to a low number of tapers.\n        ' ) #NEWLINE# N = yk . shape [ ( - 1 ) ] #NEWLINE# L = ( ( ( N / 2 ) + 1 ) if ( sides == 'onesided' ) else N ) #NEWLINE# return ( np . multiply . outer ( np . sqrt ( eigvals ) , np . ones ( L ) ) , ( 2 * K ) ) #NEWLINE# #UNINDENT# rt_eig = np . sqrt ( eigvals ) #NEWLINE# N = yk . shape [ 1 ] #NEWLINE# sdf = mtm_cross_spectrum ( yk , yk , eigvals [ : , None ] , sides = sides ) #NEWLINE# L = sdf . shape [ ( - 1 ) ] #NEWLINE# var_est = ( np . sum ( sdf , axis = ( - 1 ) ) / N ) #NEWLINE# bband_sup = ( ( 1 - eigvals ) * var_est ) #NEWLINE# sdf_iter = mtm_cross_spectrum ( yk [ : 2 ] , yk [ : 2 ] , eigvals [ : 2 , None ] , sides = sides ) #NEWLINE# err = np . zeros ( ( K , L ) ) #NEWLINE# min_pwr = ( sdf_iter . max ( ) * ( 10 ** ( ( - 150 ) / 20.0 ) ) ) #NEWLINE# default_weights = np . where ( ( sdf_iter < min_pwr ) ) [ 0 ] #NEWLINE# adaptiv_weights = np . where ( ( sdf_iter >= min_pwr ) ) [ 0 ] #NEWLINE# w_def = ( rt_eig [ : , None ] * sdf_iter [ default_weights ] ) #NEWLINE# w_def /= ( ( eigvals [ : , None ] * sdf_iter [ default_weights ] ) + bband_sup [ : , None ] ) #NEWLINE# d_sdfs = ( np . abs ( yk [ : , adaptiv_weights ] ) ** 2 ) #NEWLINE# if ( L < N ) : #NEWLINE# #INDENT# d_sdfs *= 2 #NEWLINE# #UNINDENT# sdf_iter = sdf_iter [ weights ] #NEWLINE# yk = yk [ : , adaptiv_weights ] #NEWLINE# for n in range ( max_iter ) : #NEWLINE# #INDENT# d_k = ( rt_eig [ : , None ] * sdf_iter [ None , : ] ) #NEWLINE# d_k /= ( ( eigvals [ : , None ] * sdf_iter [ None , : ] ) + bband_sup [ : , None ] ) #NEWLINE# sdf_iter = np . sum ( ( ( d_k ** 2 ) * d_sdfs ) , axis = 0 ) #NEWLINE# sdf_iter /= np . sum ( ( d_k ** 2 ) , axis = 0 ) #NEWLINE# cfn = ( eigvals [ : , None ] * ( sdf_iter [ None , : ] - d_sdfs ) ) #NEWLINE# cfn /= ( ( ( eigvals [ : , None ] * sdf_iter [ None , : ] ) + bband_sup [ : , None ] ) ** 2 ) #NEWLINE# cfn = np . sum ( cfn , axis = 0 ) #NEWLINE# if ( np . percentile ( ( cfn ** 2 ) , 95 ) < 1e-12 ) : #NEWLINE# #INDENT# break #NEWLINE# #UNINDENT# #UNINDENT# else : #NEWLINE# #INDENT# e_s = 'Breaking due to iterative meltdown in ' #NEWLINE# e_s += 'nitime.utils.adaptive_weights.' #NEWLINE# warnings . warn ( e_s , RuntimeWarning ) #NEWLINE# #UNINDENT# weights = np . zeros ( ( K , L ) ) #NEWLINE# weights [ : , adaptiv_weights ] = d_k #NEWLINE# weights [ : , default_weights ] = w_def #NEWLINE# nu = ( 2 * ( weights ** 2 ) . sum ( axis = ( - 2 ) ) ) #NEWLINE# return ( weights , nu )