<s> import warnings <EOL> import multiprocessing <EOL> from collections import OrderedDict <EOL> import numpy as np <EOL> import scipy . stats as ss <EOL> from scipy . optimize import leastsq <EOL> import period_LS_pyfftw as pLS <EOL> from feature_set import get_feature_set <EOL> from feature_set import get_feature_set_all <EOL> feature_names_list = get_feature_set ( ) <EOL> feature_names_list_all = get_feature_set_all ( ) <EOL> class ExtractFeatures : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , date , mag , err = None , n_threads = <NUM_LIT:4> , min_period = <NUM_LIT> ) : <EOL> if not isinstance ( date , np . ndarray ) : <EOL> date = np . array ( date ) <EOL> if not isinstance ( mag , np . ndarray ) : <EOL> mag = np . array ( mag ) <EOL> self . date = date <EOL> self . mag = mag <EOL> if err is not None : <EOL> if not isinstance ( err , np . ndarray ) : <EOL> err = np . array ( err ) <EOL> self . err = err <EOL> else : <EOL> self . err = np . ones ( len ( self . mag ) ) * np . std ( self . mag ) <EOL> if ( len ( self . date ) != len ( self . mag ) ) or ( len ( self . date ) != len ( self . err ) ) or ( len ( self . mag ) != len ( self . err ) ) : <EOL> raise RuntimeError ( '<STR_LIT>' ) <EOL> min_n_data = <NUM_LIT> <EOL> if len ( self . date ) < min_n_data : <EOL> warnings . warn ( '<STR_LIT>' <EOL> % min_n_data ) <EOL> n_threads = int ( n_threads ) <EOL> if n_threads > multiprocessing . cpu_count ( ) : <EOL> self . n_threads = multiprocessing . cpu_count ( ) <EOL> else : <EOL> if n_threads <= <NUM_LIT:0> : <EOL> self . n_threads = <NUM_LIT:1> <EOL> else : <EOL> self . n_threads = n_threads <EOL> min_period = float ( min_period ) <EOL> if min_period <= <NUM_LIT:0> : <EOL> self . min_period = <NUM_LIT> <EOL> else : <EOL> self . min_period = min_period <EOL> def run ( self ) : <EOL> """<STR_LIT>""" <EOL> self . shallow_run ( ) <EOL> self . deep_run ( ) <EOL> def shallow_run ( self ) : <EOL> """<STR_LIT>""" <EOL> self . n_points = len ( self . date ) <EOL> self . weight = <NUM_LIT:1.> / self . err <EOL> self . weighted_sum = np . sum ( self . weight ) <EOL> self . mean = np . mean ( self . mag ) <EOL> self . median = np . median ( self . mag ) <EOL> self . std = np . std ( self . mag ) <EOL> self . weighted_mean = np . sum ( self . mag * self . weight ) / self . weighted_sum <EOL> self . weighted_std = np . sqrt ( np . sum ( ( self . mag - self . weighted_mean ) ** <NUM_LIT:2> * self . weight ) / self . weighted_sum ) <EOL> self . skewness = ss . skew ( self . mag ) <EOL> self . kurtosis = ss . kurtosis ( self . mag ) <EOL> shapiro = ss . shapiro ( self . mag ) <EOL> self . shapiro_w = shapiro [ <NUM_LIT:0> ] <EOL> self . quartile31 = np . percentile ( self . mag , <NUM_LIT> ) - np . percentile ( self . mag , <NUM_LIT> ) <EOL> self . stetson_k = self . get_stetson_k ( self . mag , self . median , self . err ) <EOL> self . hl_amp_ratio = self . half_mag_amplitude_ratio ( <EOL> self . mag , self . median , self . weight ) <EOL> self . cusum = self . get_cusum ( self . mag ) <EOL> self . eta = self . get_eta ( self . mag , self . weighted_std ) <EOL> def deep_run ( self ) : <EOL> """<STR_LIT>""" <EOL> self . get_period_LS ( self . date , self . mag , self . n_threads , self . min_period ) <EOL> phase_folded_date = self . date % ( self . period * <NUM_LIT> ) <EOL> sorted_index = np . argsort ( phase_folded_date ) <EOL> folded_date = phase_folded_date [ sorted_index ] <EOL> folded_mag = self . mag [ sorted_index ] <EOL> self . phase_eta = self . get_eta ( folded_mag , self . weighted_std ) <EOL> self . slope_per10 , self . slope_per90 = self . slope_percentile ( folded_date , folded_mag ) <EOL> self . phase_cusum = self . get_cusum ( folded_mag ) <EOL> def get_period_LS ( self , date , mag , n_threads , min_period ) : <EOL> """<STR_LIT>""" <EOL> oversampling = <NUM_LIT> <EOL> hifac = int ( ( max ( date ) - min ( date ) ) / len ( date ) / min_period * <NUM_LIT> ) <EOL> if hifac < <NUM_LIT:100> : <EOL> hifac = <NUM_LIT:100> <EOL> fx , fy , nout , jmax , prob = pLS . fasper ( date , mag , oversampling , hifac , <EOL> n_threads ) <EOL> self . f = fx [ jmax ] <EOL> self . period = <NUM_LIT:1.> / self . f <EOL> self . period_uncertainty = self . get_period_uncertainty ( fx , fy , jmax ) <EOL> self . period_log10FAP = np . log10 ( pLS . getSignificance ( fx , fy , nout , oversampling ) [ jmax ] ) <EOL> self . period_SNR = ( fy [ jmax ] - np . median ( fy ) ) / np . std ( fy ) <EOL> order = <NUM_LIT:3> <EOL> p0 = np . ones ( order * <NUM_LIT:2> + <NUM_LIT:1> ) <EOL> date_period = ( date % self . period ) / self . period <EOL> p1 , success = leastsq ( self . residuals , p0 , <EOL> args = ( date_period , mag , order ) ) <EOL> self . amplitude = np . sqrt ( p1 [ <NUM_LIT:1> ] ** <NUM_LIT:2> + p1 [ <NUM_LIT:2> ] ** <NUM_LIT:2> ) <EOL> self . r21 = np . sqrt ( p1 [ <NUM_LIT:3> ] ** <NUM_LIT:2> + p1 [ <NUM_LIT:4> ] ** <NUM_LIT:2> ) / self . amplitude <EOL> self . r31 = np . sqrt ( p1 [ <NUM_LIT:5> ] ** <NUM_LIT:2> + p1 [ <NUM_LIT:6> ] ** <NUM_LIT:2> ) / self . amplitude <EOL> self . f_phase = np . arctan ( - p1 [ <NUM_LIT:1> ] / p1 [ <NUM_LIT:2> ] ) <EOL> self . phi21 = np . arctan ( - p1 [ <NUM_LIT:3> ] / p1 [ <NUM_LIT:4> ] ) - <NUM_LIT> * self . f_phase <EOL> self . phi31 = np . arctan ( - p1 [ <NUM_LIT:5> ] / p1 [ <NUM_LIT:6> ] ) - <NUM_LIT> * self . f_phase <EOL> """<STR_LIT>""" <EOL> def get_period_uncertainty ( self , fx , fy , jmax , fx_width = <NUM_LIT:100> ) : <EOL> """<STR_LIT>""" <EOL> start_index = jmax - fx_width <EOL> end_index = jmax + fx_width <EOL> if start_index < <NUM_LIT:0> : <EOL> start_index = <NUM_LIT:0> <EOL> if end_index > len ( fx ) - <NUM_LIT:1> : <EOL> end_index = len ( fx ) - <NUM_LIT:1> <EOL> fx_subset = fx [ start_index : end_index ] <EOL> fy_subset = fy [ start_index : end_index ] <EOL> fy_mean = np . median ( fy_subset ) <EOL> fy_std = np . std ( fy_subset ) <EOL> max_index = np . argmax ( fy_subset ) <EOL> index = np . where ( fy_subset <= fy_mean + fy_std ) [ <NUM_LIT:0> ] <EOL> left_index = index [ ( index < max_index ) ] <EOL> if len ( left_index ) == <NUM_LIT:0> : <EOL> left_index = <NUM_LIT:0> <EOL> else : <EOL> left_index = left_index [ - <NUM_LIT:1> ] <EOL> right_index = index [ ( index > max_index ) ] <EOL> if len ( right_index ) == <NUM_LIT:0> : <EOL> right_index = len ( fy_subset ) - <NUM_LIT:1> <EOL> else : <EOL> right_index = right_index [ <NUM_LIT:0> ] <EOL> half_width = ( <NUM_LIT:1.> / fx_subset [ left_index ] <EOL> - <NUM_LIT:1.> / fx_subset [ right_index ] ) / <NUM_LIT> <EOL> period_uncertainty = half_width <EOL> return period_uncertainty <EOL> def residuals ( self , pars , x , y , order ) : <EOL> """<STR_LIT>""" <EOL> return y - self . fourier_series ( pars , x , order ) <EOL> def fourier_series ( self , pars , x , order ) : <EOL> """<STR_LIT>""" <EOL> sum = pars [ <NUM_LIT:0> ] <EOL> for i in range ( order ) : <EOL> sum += pars [ i * <NUM_LIT:2> + <NUM_LIT:1> ] * np . sin ( <NUM_LIT:2> * np . pi * ( i + <NUM_LIT:1> ) * x ) + pars [ i * <NUM_LIT:2> + <NUM_LIT:2> ] * np . cos ( <NUM_LIT:2> * np . pi * ( i + <NUM_LIT:1> ) * x ) <EOL> return sum <EOL> def get_stetson_k ( self , mag , avg , err ) : <EOL> """<STR_LIT>""" <EOL> residual = ( mag - avg ) / err <EOL> stetson_k = np . sum ( np . fabs ( residual ) ) / np . sqrt ( np . sum ( residual * residual ) ) / np . sqrt ( len ( mag ) ) <EOL> return stetson_k <EOL> def half_mag_amplitude_ratio ( self , mag , avg , weight ) : <EOL> """<STR_LIT>""" <EOL> index = np . where ( mag > avg ) <EOL> lower_weight = weight [ index ] <EOL> lower_weight_sum = np . sum ( lower_weight ) <EOL> lower_mag = mag [ index ] <EOL> lower_weighted_std = np . sum ( ( lower_mag <EOL> - avg ) ** <NUM_LIT:2> * lower_weight ) / lower_weight_sum <EOL> index = np . where ( mag <= avg ) <EOL> higher_weight = weight [ index ] <EOL> higher_weight_sum = np . sum ( higher_weight ) <EOL> higher_mag = mag [ index ] <EOL> higher_weighted_std = np . sum ( ( higher_mag <EOL> - avg ) ** <NUM_LIT:2> * higher_weight ) / higher_weight_sum <EOL> return np . sqrt ( lower_weighted_std / higher_weighted_std ) <EOL> def half_mag_amplitude_ratio2 ( self , mag , avg ) : <EOL> """<STR_LIT>""" <EOL> index = np . where ( mag > avg ) <EOL> fainter_mag = mag [ index ] <EOL> lower_sum = np . sum ( ( fainter_mag - avg ) ** <NUM_LIT:2> ) / len ( fainter_mag ) <EOL> index = np . where ( mag <= avg ) <EOL> brighter_mag = mag [ index ] <EOL> higher_sum = np . sum ( ( avg - brighter_mag ) ** <NUM_LIT:2> ) / len ( brighter_mag ) <EOL> return np . sqrt ( lower_sum / higher_sum ) <EOL> def get_eta ( self , mag , std ) : <EOL> """<STR_LIT>""" <EOL> diff = mag [ <NUM_LIT:1> : ] - mag [ : len ( mag ) - <NUM_LIT:1> ] <EOL> eta = np . sum ( diff * diff ) / ( len ( mag ) - <NUM_LIT:1.> ) / std / std <EOL> return eta <EOL> def slope_percentile ( self , date , mag ) : <EOL> """<STR_LIT>""" <EOL> date_diff = date [ <NUM_LIT:1> : ] - date [ : len ( date ) - <NUM_LIT:1> ] <EOL> mag_diff = mag [ <NUM_LIT:1> : ] - mag [ : len ( mag ) - <NUM_LIT:1> ] <EOL> index = np . where ( mag_diff != <NUM_LIT:0.> ) <EOL> date_diff = date_diff [ index ] <EOL> mag_diff = mag_diff [ index ] <EOL> slope = date_diff / mag_diff <EOL> percentile_10 = np . percentile ( slope , <NUM_LIT> ) <EOL> percentile_90 = np . percentile ( slope , <NUM_LIT> ) <EOL> return percentile_10 , percentile_90 <EOL> def get_cusum ( self , mag ) : <EOL> """<STR_LIT>""" <EOL> c = np . cumsum ( mag - self . weighted_mean ) / len ( mag ) / self . weighted_std <EOL> return np . max ( c ) - np . min ( c ) <EOL> def get_features2 ( self ) : <EOL> """<STR_LIT>""" <EOL> feature_names = [ ] <EOL> feature_values = [ ] <EOL> all_vars = vars ( self ) <EOL> for name in all_vars . keys ( ) : <EOL> if not ( name == '<STR_LIT:date>' or name == '<STR_LIT>' or name == '<STR_LIT>' <EOL> or name == '<STR_LIT>' or name == '<STR_LIT>' ) : <EOL> if not ( name == '<STR_LIT:f>' or name == '<STR_LIT>' <EOL> or name == '<STR_LIT>' <EOL> or name == '<STR_LIT>' or name == '<STR_LIT>' <EOL> or name == '<STR_LIT>' or name == '<STR_LIT>' or name == '<STR_LIT>' ) : <EOL> feature_names . append ( name ) <EOL> feature_names . sort ( ) <EOL> for name in feature_names : <EOL> feature_values . append ( all_vars [ name ] ) <EOL> return feature_names , feature_values <EOL> def get_features ( self ) : <EOL> """<STR_LIT>""" <EOL> '''<STR_LIT>''' <EOL> return self . get_features_all ( ) <EOL> def get_features_all ( self ) : <EOL> """<STR_LIT>""" <EOL> features = { } <EOL> all_vars = vars ( self ) <EOL> for name in all_vars . keys ( ) : <EOL> if name in feature_names_list_all : <EOL> features [ name ] = all_vars [ name ] <EOL> features = OrderedDict ( sorted ( features . items ( ) , key = lambda t : t [ <NUM_LIT:0> ] ) ) <EOL> return features </s>
