<s> r"""<STR_LIT>""" <EOL> import numpy as np <EOL> import pandas as pd <EOL> import scipy . stats <EOL> from skbio . util . _decorator import experimental <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def closure ( mat ) : <EOL> """<STR_LIT>""" <EOL> mat = np . atleast_2d ( mat ) <EOL> if np . any ( mat < <NUM_LIT:0> ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if mat . ndim > <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if np . all ( mat == <NUM_LIT:0> , axis = <NUM_LIT:1> ) . sum ( ) > <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> mat = mat / mat . sum ( axis = <NUM_LIT:1> , keepdims = True ) <EOL> return mat . squeeze ( ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def multiplicative_replacement ( mat , delta = None ) : <EOL> r"""<STR_LIT>""" <EOL> mat = closure ( mat ) <EOL> z_mat = ( mat == <NUM_LIT:0> ) <EOL> num_feats = mat . shape [ - <NUM_LIT:1> ] <EOL> tot = z_mat . sum ( axis = - <NUM_LIT:1> , keepdims = True ) <EOL> if delta is None : <EOL> delta = ( <NUM_LIT:1.> / num_feats ) ** <NUM_LIT:2> <EOL> zcnts = <NUM_LIT:1> - tot * delta <EOL> if np . any ( zcnts ) < <NUM_LIT:0> : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> mat = np . where ( z_mat , delta , zcnts * mat ) <EOL> return mat . squeeze ( ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def perturb ( x , y ) : <EOL> r"""<STR_LIT>""" <EOL> x , y = closure ( x ) , closure ( y ) <EOL> return closure ( x * y ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def perturb_inv ( x , y ) : <EOL> r"""<STR_LIT>""" <EOL> x , y = closure ( x ) , closure ( y ) <EOL> return closure ( x / y ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def power ( x , a ) : <EOL> r"""<STR_LIT>""" <EOL> x = closure ( x ) <EOL> return closure ( x ** a ) . squeeze ( ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def inner ( x , y ) : <EOL> r"""<STR_LIT>""" <EOL> x = closure ( x ) <EOL> y = closure ( y ) <EOL> a , b = clr ( x ) , clr ( y ) <EOL> return a . dot ( b . T ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def clr ( mat ) : <EOL> r"""<STR_LIT>""" <EOL> mat = closure ( mat ) <EOL> lmat = np . log ( mat ) <EOL> gm = lmat . mean ( axis = - <NUM_LIT:1> , keepdims = True ) <EOL> return ( lmat - gm ) . squeeze ( ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def clr_inv ( mat ) : <EOL> r"""<STR_LIT>""" <EOL> return closure ( np . exp ( mat ) ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def ilr ( mat , basis = None , check = True ) : <EOL> r"""<STR_LIT>""" <EOL> mat = closure ( mat ) <EOL> if basis is None : <EOL> basis = clr_inv ( _gram_schmidt_basis ( mat . shape [ - <NUM_LIT:1> ] ) ) <EOL> elif check : <EOL> _check_orthogonality ( basis ) <EOL> return inner ( mat , basis ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def ilr_inv ( mat , basis = None , check = True ) : <EOL> r"""<STR_LIT>""" <EOL> if basis is None : <EOL> basis = _gram_schmidt_basis ( mat . shape [ - <NUM_LIT:1> ] + <NUM_LIT:1> ) <EOL> elif check : <EOL> _check_orthogonality ( basis ) <EOL> return clr_inv ( np . dot ( mat , basis ) ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def centralize ( mat ) : <EOL> r"""<STR_LIT>""" <EOL> mat = closure ( mat ) <EOL> cen = scipy . stats . gmean ( mat , axis = <NUM_LIT:0> ) <EOL> return perturb_inv ( mat , cen ) <EOL> @ experimental ( as_of = "<STR_LIT>" ) <EOL> def ancom ( table , grouping , <EOL> alpha = <NUM_LIT> , <EOL> tau = <NUM_LIT> , <EOL> theta = <NUM_LIT:0.1> , <EOL> multiple_comparisons_correction = None , <EOL> significance_test = None ) : <EOL> r"""<STR_LIT>""" <EOL> if not isinstance ( table , pd . DataFrame ) : <EOL> raise TypeError ( '<STR_LIT>' <EOL> '<STR_LIT>' % type ( table ) . __name__ ) <EOL> if not isinstance ( grouping , pd . Series ) : <EOL> raise TypeError ( '<STR_LIT>' <EOL> '<STR_LIT>' % type ( grouping ) . __name__ ) <EOL> if np . any ( table <= <NUM_LIT:0> ) : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' <EOL> ) <EOL> if not <NUM_LIT:0> < alpha < <NUM_LIT:1> : <EOL> raise ValueError ( '<STR_LIT>' % alpha ) <EOL> if not <NUM_LIT:0> < tau < <NUM_LIT:1> : <EOL> raise ValueError ( '<STR_LIT>' % tau ) <EOL> if not <NUM_LIT:0> < theta < <NUM_LIT:1> : <EOL> raise ValueError ( '<STR_LIT>' % theta ) <EOL> if multiple_comparisons_correction is not None : <EOL> if multiple_comparisons_correction != '<STR_LIT>' : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' <EOL> % multiple_comparisons_correction ) <EOL> if ( grouping . isnull ( ) ) . any ( ) : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> if ( table . isnull ( ) ) . any ( ) . any ( ) : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> groups , _grouping = np . unique ( grouping , return_inverse = True ) <EOL> grouping = pd . Series ( _grouping , index = grouping . index ) <EOL> num_groups = len ( groups ) <EOL> if num_groups == len ( grouping ) : <EOL> raise ValueError ( <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> if num_groups == <NUM_LIT:1> : <EOL> raise ValueError ( <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" <EOL> "<STR_LIT>" ) <EOL> if significance_test is None : <EOL> significance_test = scipy . stats . f_oneway <EOL> table_index_len = len ( table . index ) <EOL> grouping_index_len = len ( grouping . index ) <EOL> mat , cats = table . align ( grouping , axis = <NUM_LIT:0> , join = '<STR_LIT>' ) <EOL> if ( len ( mat ) != table_index_len or len ( cats ) != grouping_index_len ) : <EOL> raise ValueError ( '<STR_LIT>' <EOL> '<STR_LIT>' ) <EOL> n_feat = mat . shape [ <NUM_LIT:1> ] <EOL> _logratio_mat = _log_compare ( mat . values , cats . values , significance_test ) <EOL> logratio_mat = _logratio_mat + _logratio_mat . T <EOL> if multiple_comparisons_correction == '<STR_LIT>' : <EOL> logratio_mat = np . apply_along_axis ( _holm_bonferroni , <EOL> <NUM_LIT:1> , logratio_mat ) <EOL> np . fill_diagonal ( logratio_mat , <NUM_LIT:1> ) <EOL> W = ( logratio_mat < alpha ) . sum ( axis = <NUM_LIT:1> ) <EOL> c_start = W . max ( ) / n_feat <EOL> if c_start < theta : <EOL> reject = np . zeros_like ( W , dtype = bool ) <EOL> else : <EOL> cutoff = c_start - np . linspace ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT:5> ) <EOL> prop_cut = np . array ( [ ( W > n_feat * cut ) . mean ( ) for cut in cutoff ] ) <EOL> dels = np . abs ( prop_cut - np . roll ( prop_cut , - <NUM_LIT:1> ) ) <EOL> dels [ - <NUM_LIT:1> ] = <NUM_LIT:0> <EOL> if ( dels [ <NUM_LIT:0> ] < tau ) and ( dels [ <NUM_LIT:1> ] < tau ) and ( dels [ <NUM_LIT:2> ] < tau ) : <EOL> nu = cutoff [ <NUM_LIT:1> ] <EOL> elif ( dels [ <NUM_LIT:0> ] >= tau ) and ( dels [ <NUM_LIT:1> ] < tau ) and ( dels [ <NUM_LIT:2> ] < tau ) : <EOL> nu = cutoff [ <NUM_LIT:2> ] <EOL> elif ( dels [ <NUM_LIT:1> ] >= tau ) and ( dels [ <NUM_LIT:2> ] < tau ) and ( dels [ <NUM_LIT:3> ] < tau ) : <EOL> nu = cutoff [ <NUM_LIT:3> ] <EOL> else : <EOL> nu = cutoff [ <NUM_LIT:4> ] <EOL> reject = ( W >= nu * n_feat ) <EOL> labs = mat . columns <EOL> return pd . DataFrame ( { '<STR_LIT>' : pd . Series ( W , index = labs ) , <EOL> '<STR_LIT>' : pd . Series ( reject , index = labs ) } ) <EOL> def _holm_bonferroni ( p ) : <EOL> """<STR_LIT>""" <EOL> K = len ( p ) <EOL> sort_index = - np . ones ( K , dtype = np . int64 ) <EOL> sorted_p = np . sort ( p ) <EOL> sorted_p_adj = sorted_p * ( K - np . arange ( K ) ) <EOL> for j in range ( K ) : <EOL> idx = ( p == sorted_p [ j ] ) & ( sort_index < <NUM_LIT:0> ) <EOL> num_ties = len ( sort_index [ idx ] ) <EOL> sort_index [ idx ] = np . arange ( j , ( j + num_ties ) , dtype = np . int64 ) <EOL> sorted_holm_p = [ min ( [ max ( sorted_p_adj [ : k ] ) , <NUM_LIT:1> ] ) <EOL> for k in range ( <NUM_LIT:1> , K + <NUM_LIT:1> ) ] <EOL> holm_p = [ sorted_holm_p [ sort_index [ k ] ] for k in range ( K ) ] <EOL> return holm_p <EOL> def _log_compare ( mat , cats , <EOL> significance_test = scipy . stats . ttest_ind ) : <EOL> """<STR_LIT>""" <EOL> r , c = mat . shape <EOL> log_ratio = np . zeros ( ( c , c ) ) <EOL> log_mat = np . log ( mat ) <EOL> cs = np . unique ( cats ) <EOL> def func ( x ) : <EOL> return significance_test ( * [ x [ cats == k ] for k in cs ] ) <EOL> for i in range ( c - <NUM_LIT:1> ) : <EOL> ratio = ( log_mat [ : , i ] . T - log_mat [ : , i + <NUM_LIT:1> : ] . T ) . T <EOL> m , p = np . apply_along_axis ( func , <EOL> axis = <NUM_LIT:0> , <EOL> arr = ratio ) <EOL> log_ratio [ i , i + <NUM_LIT:1> : ] = np . squeeze ( np . array ( p . T ) ) <EOL> return log_ratio <EOL> def _gram_schmidt_basis ( n ) : <EOL> """<STR_LIT>""" <EOL> basis = np . zeros ( ( n , n - <NUM_LIT:1> ) ) <EOL> for j in range ( n - <NUM_LIT:1> ) : <EOL> i = j + <NUM_LIT:1> <EOL> e = np . array ( [ ( <NUM_LIT:1> / i ) ] * i + [ - <NUM_LIT:1> ] + <EOL> [ <NUM_LIT:0> ] * ( n - i - <NUM_LIT:1> ) ) * np . sqrt ( i / ( i + <NUM_LIT:1> ) ) <EOL> basis [ : , j ] = e <EOL> return basis . T <EOL> def _check_orthogonality ( basis ) : <EOL> """<STR_LIT>""" <EOL> if not np . allclose ( inner ( basis , basis ) , np . identity ( len ( basis ) ) , <EOL> rtol = <NUM_LIT> , atol = <NUM_LIT> ) : <EOL> raise ValueError ( "<STR_LIT>" ) </s>
