<s> """<STR_LIT>""" <EOL> from numpy import * <EOL> try : <EOL> import pyfftw <EOL> is_pyfftw = True <EOL> except : <EOL> from numpy . fft import * <EOL> is_pyfftw = False <EOL> def __spread__ ( y , yy , n , x , m ) : <EOL> """<STR_LIT>""" <EOL> nfac = [ <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:2> , <NUM_LIT:6> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ] <EOL> if m > <NUM_LIT> : <EOL> print '<STR_LIT>' <EOL> return <EOL> ix = long ( x ) <EOL> if x == float ( ix ) : <EOL> yy [ ix ] = yy [ ix ] + y <EOL> else : <EOL> ilo = long ( x - <NUM_LIT:0.5> * float ( m ) + <NUM_LIT:1.0> ) <EOL> ilo = min ( max ( ilo , <NUM_LIT:1> ) , n - m + <NUM_LIT:1> ) <EOL> ihi = ilo + m - <NUM_LIT:1> <EOL> nden = nfac [ m ] <EOL> fac = x - ilo <EOL> for j in range ( ilo + <NUM_LIT:1> , ihi + <NUM_LIT:1> ) : fac = fac * ( x - j ) <EOL> yy [ ihi ] = yy [ ihi ] + y * fac / ( nden * ( x - ihi ) ) <EOL> for j in range ( ihi - <NUM_LIT:1> , ilo - <NUM_LIT:1> , - <NUM_LIT:1> ) : <EOL> nden = ( nden / ( j + <NUM_LIT:1> - ilo ) ) * ( j - ihi ) <EOL> yy [ j ] = yy [ j ] + y * fac / ( nden * ( x - j ) ) <EOL> def fasper ( x , y , ofac , hifac , n_threads , MACC = <NUM_LIT:4> ) : <EOL> """<STR_LIT>""" <EOL> n = long ( len ( x ) ) <EOL> if n != len ( y ) : <EOL> print '<STR_LIT>' <EOL> return <EOL> nout = <NUM_LIT:0.5> * ofac * hifac * n <EOL> nfreqt = long ( ofac * hifac * n * MACC ) <EOL> nfreq = <NUM_LIT:64> L <EOL> while nfreq < nfreqt : <EOL> nfreq = <NUM_LIT:2> * nfreq <EOL> ndim = long ( <NUM_LIT:2> * nfreq ) <EOL> ave = y . mean ( ) <EOL> var = ( ( y - y . mean ( ) ) ** <NUM_LIT:2> ) . sum ( ) / ( len ( y ) - <NUM_LIT:1> ) <EOL> xmin = x . min ( ) <EOL> xmax = x . max ( ) <EOL> xdif = xmax - xmin <EOL> if is_pyfftw : <EOL> wk1 = pyfftw . n_byte_align_empty ( int ( ndim ) , <NUM_LIT:16> , '<STR_LIT>' ) * <NUM_LIT:0.> <EOL> wk2 = pyfftw . n_byte_align_empty ( int ( ndim ) , <NUM_LIT:16> , '<STR_LIT>' ) * <NUM_LIT:0.> <EOL> else : <EOL> wk1 = zeros ( ndim , dtype = '<STR_LIT>' ) <EOL> wk2 = zeros ( ndim , dtype = '<STR_LIT>' ) <EOL> fac = ndim / ( xdif * ofac ) <EOL> fndim = ndim <EOL> ck = ( ( x - xmin ) * fac ) % fndim <EOL> ckk = ( <NUM_LIT> * ck ) % fndim <EOL> for j in range ( <NUM_LIT:0> L , n ) : <EOL> __spread__ ( y [ j ] - ave , wk1 , ndim , ck [ j ] , MACC ) <EOL> __spread__ ( <NUM_LIT:1.0> , wk2 , ndim , ckk [ j ] , MACC ) <EOL> if is_pyfftw : <EOL> fft_wk1 = pyfftw . builders . ifft ( wk1 , planner_effort = '<STR_LIT>' , <EOL> threads = n_threads ) <EOL> wk1 = fft_wk1 ( ) * len ( wk1 ) <EOL> fft_wk2 = pyfftw . builders . ifft ( wk2 , planner_effort = '<STR_LIT>' , <EOL> threads = n_threads ) <EOL> wk2 = fft_wk2 ( ) * len ( wk2 ) <EOL> else : <EOL> wk1 = ifft ( wk1 ) * len ( wk1 ) <EOL> wk2 = ifft ( wk2 ) * len ( wk1 ) <EOL> wk1 = wk1 [ <NUM_LIT:1> : nout + <NUM_LIT:1> ] <EOL> wk2 = wk2 [ <NUM_LIT:1> : nout + <NUM_LIT:1> ] <EOL> rwk1 = wk1 . real <EOL> iwk1 = wk1 . imag <EOL> rwk2 = wk2 . real <EOL> iwk2 = wk2 . imag <EOL> df = <NUM_LIT:1.0> / ( xdif * ofac ) <EOL> hypo2 = <NUM_LIT> * abs ( wk2 ) <EOL> hc2wt = rwk2 / hypo2 <EOL> hs2wt = iwk2 / hypo2 <EOL> cwt = sqrt ( <NUM_LIT:0.5> + hc2wt ) <EOL> swt = sign ( hs2wt ) * ( sqrt ( <NUM_LIT:0.5> - hc2wt ) ) <EOL> den = <NUM_LIT:0.5> * n + hc2wt * rwk2 + hs2wt * iwk2 <EOL> cterm = ( cwt * rwk1 + swt * iwk1 ) ** <NUM_LIT> / den <EOL> sterm = ( cwt * iwk1 - swt * rwk1 ) ** <NUM_LIT> / ( n - den ) <EOL> wk1 = df * ( arange ( nout , dtype = '<STR_LIT:float>' ) + <NUM_LIT:1.> ) <EOL> wk2 = ( cterm + sterm ) / ( <NUM_LIT> * var ) <EOL> pmax = wk2 . max ( ) <EOL> jmax = wk2 . argmax ( ) <EOL> expy = exp ( - pmax ) <EOL> effm = <NUM_LIT> * ( nout ) / ofac <EOL> prob = effm * expy <EOL> if prob > <NUM_LIT> : <EOL> prob = <NUM_LIT:1.0> - ( <NUM_LIT:1.0> - expy ) ** effm <EOL> return wk1 , wk2 , nout , jmax , prob <EOL> def getSignificance ( wk1 , wk2 , nout , ofac ) : <EOL> """<STR_LIT>""" <EOL> expy = exp ( - wk2 ) <EOL> effm = <NUM_LIT> * ( nout ) / ofac <EOL> sig = effm * expy <EOL> ind = ( sig > <NUM_LIT> ) . nonzero ( ) <EOL> sig [ ind ] = <NUM_LIT:1.0> - ( <NUM_LIT:1.0> - expy [ ind ] ) ** effm <EOL> return sig </s>
