# Dependencies import os import sys import subprocess subprocess.check_call('pip install -r /opt/ml/processing/input/dependencies/requirements-2.txt', shell=True) import logging logging.basicConfig( level=logging.INFO, format='%(asctime)s.%(msecs)03d %(levelname)s %(module)s - %(funcName)s: %(message)s', datefmt='%Y-%m-%d %H:%M:%S', ) import pandas as pd import boto3 import math import random import scipy import numpy as np import statsmodels.formula.api as smf import statsmodels.api as sm import pmdarima as pm import time from datetime import datetime, date, timedelta from statsmodels.tsa.arima.model import ARIMA from statsmodels.tools.eval_measures import mse, rmse from statsmodels.tsa.statespace.tools import diff from scipy import fftpack from multiprocessing import Pool, cpu_count def timing_val(func): ''' Utility function that calculates the execution time of other functions. When @timing_val is present, the function will return a tuple containing result and execution time. From source: http://www.daniweb.com/code/snippet368.html ''' def wrapper(*arg, **kw): t1 = time.time() res = func(*arg, **kw) t2 = time.time() return (t2 - t1), res, func.__name__ return wrapper ################################################################################################################## ############################################# FOURIER ############################################################ ################################################################################################################## def fourier_table(one_track_df, loop=60): ''' Apply Fast Fourier Transform (FFT) to individual time series. FFT is a method for expressing a function as a sum of its periodic components. This is used to clean the noise from the data using low pass filter, which cuts off high frequencies. The function also identifies the inflection points of the FFT time series. one_track_df -- DataFrame containing the time series data loop -- the number of loops for the Fourier Transform. ''' # Computing FFT and applying low pass filter y = np.array(one_track_df['STREAMS']) y_fft_filtered = fftpack.fft(y).copy() freq = fftpack.fftfreq(len(y), d=1 / len(y)) # Filtering the frequencies cut_off = len(y) / loop y_fft_filtered[np.abs(freq) > cut_off] = 0 Fourier = fftpack.ifft(y_fft_filtered) # Computing second derivative to identify inflection points y_diff = pd.Series(Fourier).diff() y_diff_diff = y_diff.diff() locs2 = abs(np.diff(np.sign(y_diff_diff))) zero_diff_diff = np.where(np.logical_and(locs2 != 0, np.isnan(locs2) == False)) inflection_point = one_track_df.iloc[zero_diff_diff[0]] # Compiling table containing FFT columns Fourier = pd.DataFrame(Fourier, columns=['Fourier']) Fourier['ACTIVITY_DATE'] = np.array(one_track_df['ACTIVITY_DATE']) streams = pd.DataFrame(one_track_df[['STREAMS', 'ACTIVITY_DATE']]) df_fourier = streams.merge(Fourier, on=['ACTIVITY_DATE']) df_fourier['Inflection_Point'] = np.where(df_fourier['ACTIVITY_DATE'] .isin(inflection_point['ACTIVITY_DATE']), 1, 0) df_fourier['Fourier_real_part'] = np.array(df_fourier['Fourier']).real df_fourier['ISRC_KEY'] = np.array(one_track_df['ISRC_KEY']) df_fourier['ISRC'] = np.array(one_track_df['ISRC']) df_fourier['TRANSACTION_COUNTRY_CODE'] = np.array(one_track_df['TRANSACTION_COUNTRY_CODE']) return df_fourier @timing_val def compile_fourier_table(df): '''Assembling the dataframe containing FFT and Inflection Point columns''' res = pd.DataFrame(columns=['ACTIVITY_DATE', 'STREAMS', 'Fourier', 'Inflection_Point', 'ISRC_KEY', 'ISRC', 'TRANSACTION_COUNTRY_CODE']) unique_isrcs = df['ISRC_KEY'].unique() for isrc in unique_isrcs: one_track_df = unique_isrc_df(isrc, df) if check_one_track_df(one_track_df): one_track_fourier = fourier_table(one_track_df, loop=60) res = pd.concat([res, one_track_fourier]) else: logging.error(f"Error on {isrc}: data has length {len(one_track_df)}") return res #def compute_fourier_for_chunk(chunk_df, folder, run_id, chunk_id, dry_run): def compute_fourier_for_chunk(chunk_df, chunk_id, par_id): ''' Compute FastFourierTransform for a the given chunk of data, and save the data to S3 chunk_df -- current dataframe for running chunk folder -- where work is stored run_id -- ID of current run chunk_id -- ID of current chunk of data running par_id -- parallel job number dry_run -- bolean flag indicating whether the results are saved. ''' logging.info(f"Starting to compute FFT on data chunk {chunk_id} and parallel job {par_id,} ...") fourier_table = compile_fourier_table(chunk_df) logging.info(f"FFT on data chunk {chunk_id} and parallel job {par_id,} finished.") fourier_df = pd.DataFrame(fourier_table[1]) fourier_df['ACTIVITY_DATE'] = pd.to_datetime(fourier_df['ACTIVITY_DATE']).dt.date logging.info(f"Saving FFT data to S3 ...") # save_fourier_s3(fourier_df, folder, 'fourier_table', run_id, chunk_id, dry_run) # logging.info(f"FFT data saved.") return fourier_df def save_fourier_s3(df, folder, prefix, run_id, chunk_id, par_id, dry_run): ''' Saving FFT calculated data to S3 df -- dataframe to save prefix -- assigned prefix for file name folder -- where work is stored run_id -- ID of current run chunk_id -- ID of current chunk of data running dry_run -- bolean flag indicating whether the results are saved. ''' if dry_run: logging.warning(f"DRY_RUN flag is True, no FFT data will be saved.") else: s3 = boto3.client('s3') bucket_name = 'dev-cucumbers' filepath = "{}/{}_run_{}_chunk_{}_par{}.csv".format(folder, prefix, run_id, chunk_id, par_id) csv_buffer = df.to_csv(index=False).encode('utf-8') s3.put_object(Body=csv_buffer, Bucket=bucket_name, Key=filepath) logging.info(f"Table saved to S3 bucket: {bucket_name}, with file name: {filepath}") def save_fourier_s3_inflection_only(df, folder, prefix, run_id, chunk_id, par_id, dry_run): if dry_run: logging.warning(f"DRY_RUN flag is True, no FFT data will be saved.") else: s3 = boto3.client('s3') bucket_name = 'dev-cucumbers' filepath = "{}/{}_run_{}_chunk_{}_par{}.csv".format(folder, prefix, run_id, chunk_id, par_id) csv_buffer = df.to_csv(index=False).encode('utf-8') s3.put_object(Body=csv_buffer, Bucket=bucket_name, Key=filepath) logging.info(f"Table saved to S3 bucket: {bucket_name}, with file name: {filepath}") ################################################################################################################## ############################################# ARIMA ############################################################## ################################################################################################################## @timing_val def take_isrcs_with_inflection_last_week(fourier_df): ''' Function that selects only ISRCs that recorded an inflection point in the last 7 days of their lifecycle ''' first_day_pred = fourier_df['ACTIVITY_DATE'].max() - timedelta(days=7) list_for_pred = fourier_df[(fourier_df['Inflection_Point'] == 1) & (fourier_df['ACTIVITY_DATE'] > first_day_pred)]['ISRC_KEY'].unique().tolist() subset_for_pred = fourier_df[fourier_df['ISRC_KEY'].isin(list_for_pred)].copy() return subset_for_pred @timing_val def take_isrcs_with_inflection_last_week_GB(fourier_df): ''' Function that selects only ISRCs that recorded an inflection point in the last 7 days of their lifecycle ''' first_day_pred = fourier_df['ACTIVITY_DATE'].max() - timedelta(days=7) list_for_pred = fourier_df[(fourier_df['Inflection_Point'] == 1) & (fourier_df['ACTIVITY_DATE'] > first_day_pred)]['ISRC_KEY'].unique().tolist() subset_for_pred = fourier_df[fourier_df['ISRC_KEY'].isin(list_for_pred)].copy() subset_for_pred = subset_for_pred[subset_for_pred['TRANSACTION_COUNTRY_CODE'] == 'GB'].copy() return subset_for_pred def iterate_group_by_key(df, col_for_key, sort_data=True): ''' Optimized function (generator) that iterates over the DataFrame and groups it by a specified column. It operates on a sorted dataframe, so all rows for a given key (here ISRC) are together. It yields one group at a time as a new DataFrame. ''' if sort_data: df = df.sort_values(by=col_for_key) # helper function that returns the value in the column specified by col_for_key at the index i. def key_at(i): return np.array(df[col_for_key].iloc[i]) index, size = (0, df.shape[0]) while index < size: current_key = key_at(index) res = [] while index < size and list(current_key) == list(key_at(index)): res.append(df.iloc[index]) index = index + 1 resdf = pd.DataFrame(res, columns=df.columns) yield resdf @timing_val def auto_arima(df_isrc): ''' ARIMA model that applies a stepwise selection to decide autoregressive parameters. Autoregressive parameters: e.g ARIMA(p, d, q) - p refers to an autoregressive process AR(p): Y_t = α_1Y_(t_1) + ... + α_pY_(t_p) + Z_t, where α_1 is the coefficient of the AR process and Y_t regresses in its past values Y_(t_1), and Z_t is the purely random process (error term). Effectively, 'p' stands for the number of lagged terms used in the autoregressive component. - d refers to the d_th order differenced process: integration makes the process stationary when non stationary. When d=0, process is stationary (mean and variance don't change over time). Effectively, 'd' indicates how many times the time series is differenced before reaching stationarity. - q refers to a moving average process MA(q): Y_t = λ_0Z_t + λ_1Z_(t-1) +... + λ_qZ_(t-q), where Z_t (error term) is the purely random process and is given a coefficient of 1 (λ_0=1). Effectively, 'q' indicates the number of lagged error terms used in the moving average component. For further documentation on pmdarima.auto_arima() see: https://alkaline-ml.com/pmdarima/modules/generated/pmdarima.arima.auto_arima.html ''' cutoff_train_test_sets = 133 (train, test) = (df_isrc.iloc[:cutoff_train_test_sets], df_isrc.iloc[cutoff_train_test_sets:]) train.index = pd.to_datetime(train['ACTIVITY_DATE']) train = train.sort_index(axis=0) auto_df = train[['STREAMS']].copy() auto_model = pm.auto_arima(auto_df, start_p=1, start_q=1, test='adf', max_p=3, max_q=3, m=7, start_P=0, seasonal=True, d=None, D=1, trace=False, error_action='ignore', trend='ct', suppress_warnings=True, # stationary=False, stepwise=True) return (auto_model.seasonal_order, auto_model.order) @timing_val def arima_iscrs(df): ''' Function that applies the winning model from auto_arima() and applies cross validation: it trains on the first 133 days of the series (train set) and forecasts the last 7 days (test set). The function returns a table with model evaluation values: - pred_type: whether actual streams are above/below/crossing predicted streams - mse (mean squared error): average squared difference between observed and predicted values (measures of error in statistical models) - avg_streams_train: mean of total streams in train set - avg_streams_test: mean of total streams in test set - median_streams_train: median of total streams in train set - median_streams_test: median of total streams in test set - linear_gradient_train: slope of train set - linear_gradient_test: slope of test set - sum_forecast_errors: sum of all forecast errors in test set. A forecast error is the difference between an observed value and its predicted value. - len_df: length of the time series. ''' # Setting up empty dataframe new_cs_df = pd.DataFrame(columns=['ISRC_KEY', 'ISRC', 'pred_type', 'mse', 'avg_streams_train', 'avg_streams_test', 'median_streams_train', 'median_streams_test', 'linear_gradient_train', 'linear_gradient_test', 'sum_forecast_errors', 'len_df', 'earliest_inf_point', 'TRANSACTION_COUNTRY_CODE']) new_ts_df = pd.DataFrame(columns=['ISRC_KEY', 'ISRC', 'ACTIVITY_DATE', 'STREAMS', 'predicted_values', 'len_df', 'TRANSACTION_COUNTRY_CODE']) count = 1 total_isrcs = df['ISRC_KEY'].nunique() auto_arima_durations = [] # Subsetting by individual ISRC dataframe (for univariate analysis) for subdf in iterate_group_by_key(df, ['ISRC_KEY']): df_isrc = subdf isrc_key = df_isrc['ISRC_KEY'].iloc[0] isrc = df_isrc['ISRC'].iloc[0] df_isrc = df_isrc.sort_values(by=['ACTIVITY_DATE']) activity = df_isrc['ACTIVITY_DATE'] streams = df_isrc['STREAMS'] country = df_isrc['TRANSACTION_COUNTRY_CODE'].iloc[0] latest_date = df_isrc['ACTIVITY_DATE'].max() date_7_days_ago = latest_date - pd.to_timedelta(6, unit='d') df_last_7_days = df_isrc[df_isrc['ACTIVITY_DATE'] >= date_7_days_ago] filtered_df = df_last_7_days[df_last_7_days['Inflection_Point'] == 1] first_inflection_point = filtered_df['ACTIVITY_DATE'].min() try: # Setting up cross-validation cutoff_train_test_sets = 133 start_pred = 133 end_pred = 140 (train, test) = (df_isrc.iloc[:cutoff_train_test_sets], df_isrc.iloc[cutoff_train_test_sets:]) (end_test, end_train) = (len(test), len(train)) # print(len(test)) # print(len(train)) # Modelling happens here _, (seasonal_order, order), _ = auto_arima(train) model = ARIMA(train['STREAMS'], order=order, seasonal_order=seasonal_order, enforce_stationarity=False).fit() pred = model.predict(start=end_train, end=(end_train + end_test - 1)) # pred = model.predict(n_periods=7, return_conf_int=False) # in_sample_pred = model.get_prediction(start= 1, end = (end_train + end_test-1), dynamic=False) # predicted_values = in_sample_pred.predicted_mean[:end_pred] predicted_values = model.predict(start=0, end=139) # print('pred 7 days: ', pred) # # print(in_sample_pred) # print('pred 140 days: ', predicted_values) # print('streams: ',test['STREAMS']) # Extracting relevant metrics forecast_errors = np.subtract(np.array(test['STREAMS']), pred) # print('FE: ',forecast_errors) mse = np.square(forecast_errors).mean() avg_streams_train = train['STREAMS'].mean() avg_streams_test = test['STREAMS'].mean() median_streams_train = train['STREAMS'].median() median_streams_test = test['STREAMS'].median() linear_gradient_train = (train.iloc[-1]['STREAMS'] - train.iloc[0, test.columns.get_loc('STREAMS')]) / len( train) linear_gradient_test = (test.iloc[-1]['STREAMS'] - test.iloc[0, test.columns.get_loc('STREAMS')]) / len( test) sum_forecast_errors = round(forecast_errors.sum(), 3) all_positives = all(map(lambda x: x > 0, forecast_errors)) all_negatives = all(map(lambda x: x < 0, forecast_errors)) pred_categorical = None if all_positives: pred_categorical = 'actuals_above_predicted' elif all_negatives: pred_categorical = 'actuals_below_predicted' else: pred_categorical = 'actuals_crossing_predicted' # Populating cross-sectional dataframe new_cs_df = pd.concat([new_cs_df, pd.DataFrame([{ 'ISRC': isrc, 'ISRC_KEY': isrc_key, 'pred_type': pred_categorical, 'mse': mse, 'avg_streams_train': avg_streams_train, 'avg_streams_test': avg_streams_test, 'median_streams_train': median_streams_train, 'median_streams_test': median_streams_test, 'linear_gradient_train': linear_gradient_train, 'linear_gradient_test': linear_gradient_test, 'sum_forecast_errors': sum_forecast_errors, 'len_df': df_isrc.shape[0], 'earliest_inf_point': first_inflection_point, 'TRANSACTION_COUNTRY_CODE': country }], columns=new_cs_df.columns)]) # Populating time-series dataframe isrcs_col = np.full( shape=df_isrc.shape[0], fill_value=isrc_key, dtype=object ) len_df_col = np.full( shape=df_isrc.shape[0], fill_value=df_isrc.shape[0], dtype=int ) new_ts_df = pd.concat([new_ts_df, pd.DataFrame({ 'ISRC_KEY': isrcs_col, 'ISRC': isrc, 'ACTIVITY_DATE': activity, 'STREAMS': streams, 'predicted_values': predicted_values, 'len_df': len_df_col, 'TRANSACTION_COUNTRY_CODE': country }, columns=new_ts_df.columns)]) if count % 100 == 0: logging.info(f"Parallel worker processed {count} out of {total_isrcs} ISRCs.") count = count + 1 except Exception as e: logging.error(f"Error while processing ISRC {isrc}: '{e}'") # Computing avg calculation time for auto-arima (and make sure we don't divide by 0 !) avg_duration = sum(auto_arima_durations) / max(1, len(auto_arima_durations)) logging.info(f"Average auto-arima duration was {avg_duration} seconds for {total_isrcs} ISRCs.") #logging.info(f"ISRCs: {list(isrc)}") return (new_cs_df, new_ts_df) def save_chunk_s3(df, folder, prefix, run_id, chunk_id, par_id, dry_run): ''' Save the given parallel chunk result into S3 df -- dataframe to save prefix -- assigned prefix for file name folder -- where work is stored run_id -- ID of current run chunk_id -- ID of current chunk of data running dry_run -- bolean flag indicating whether the results are saved. ''' if dry_run: logging.warning(f"DRY_RUN flag is True, no parallel chunk will be saved.") else: s3 = boto3.client('s3') bucket_name = 'dev-cucumbers' filepath = "{}/{}_run_{}_chunk_{}_par_{}.csv".format(folder, prefix, run_id, chunk_id, par_id) csv_buffer = df.to_csv(index=False).encode('utf-8') s3.put_object(Body=csv_buffer, Bucket=bucket_name, Key=filepath) logging.info(f"Table saved to S3 bucket: {bucket_name}, with file name: {filepath}") ################################################################################################################## ############################################# UTILS ############################################################## ################################################################################################################## def check_one_track_df(one_track_df): '''Checking that the length of the individual time series is 140 days''' return len(one_track_df) == 140 def unique_isrc_df(isrc, df): ''' As time series analysis is univariate, decomposing the dataframe as to have one-track dataframes at time of analysis. ''' one_track_df = df[df['ISRC_KEY'] == isrc].copy() one_track_df['ACTIVITY_DATE'] = pd.to_datetime(one_track_df['ACTIVITY_DATE']) one_track_df['ACTIVITY_DATE'] = one_track_df['ACTIVITY_DATE'].dt.date one_track_df = one_track_df[['ISRC_KEY', 'ISRC', 'ACTIVITY_DATE', 'STREAMS', 'TRANSACTION_COUNTRY_CODE']] return one_track_df.sort_values(by='ACTIVITY_DATE') def data_for_chunk(df, chunk_id, chunk_count): ''' Function sorting ISRCs and splitting them into n chuncks and taking only one of them (chunk_id) ''' IDs = df['ISRC_KEY'].unique().tolist() IDs.sort() chunk_size = len(IDs) // chunk_count chunked = [IDs[n:n + chunk_size] for n in range(0, len(IDs), chunk_size)] ids = chunked[chunk_id] chunk_df = df[df['ISRC_KEY'].isin(ids)] return chunk_df def split_universe(df, n): ''' Function splitting dataframe for parallelism df -- initial dataframe n -- number of chunks ''' # list of sorted ISRCs into n chunks def chunks(l, n): for i in range(0, len(l), n): yield l[i:i + n] isrcs = df['ISRC_KEY'].unique().tolist() isrcs.sort() isrc_chunks = chunks(isrcs, (len(isrcs) // n) + 1) return list(isrc_chunks) @timing_val def split_work(full_df, n, f, folder, run_id, chunk_id, dry_run): ''' full_df -- full dataframe that will be split into n parallel tasks n -- number of parallel tasks f -- function handling one parallel task folder -- where work is stored run_id -- ID of current run chunk_id -- ID of current chunk of data running dry_run -- bolean flag indicating whether the results are saved. ''' isrc_chunks = split_universe(full_df, n) splitted_df = [] # For each chunk, creating an individual dataframe that only includes rows with the same ID for chunk in isrc_chunks: df_chunk = full_df[full_df['ISRC_KEY'].isin(chunk)] splitted_df.append(df_chunk) # Unique run ID generator based on timestamp logging.info("ISRCs have been splitted in {} dataframes for run id={}".format(n, run_id)) # Using multiprocessing pool to apply the function to each chunk with Pool(n) as p: # Adding indices to the list of chunks splitted_df_with_index = list(enumerate(splitted_df)) # List of arguments for the function splitted_df_with_args = map(lambda x: (dry_run, folder, run_id, chunk_id, x[0], x[1]), splitted_df_with_index) p.map(f, splitted_df_with_args) logging.info("Moments calculation finished.") def work_load(arg): ''' Function that returns the result of arima_iscrs() and saves it in S3''' (dry_run, folder, run_id, chunk_id, par_id, df) = arg logging.info( f"Arguments folder={folder}, run-id={run_id}, chunk-id={chunk_id}, par-id={par_id}, df type={type(df)}") # Fourier logging.info(f"Calculation chunk {CHUNK_ID} subset contains {df.shape[0]} rows.") logging.info(f"Splitting work into {parallelism} parts for parallel processing") # Compute FastFourierTransform and save the data to S3 fourier_df = compute_fourier_for_chunk(df, chunk_id, par_id) save_fourier_s3(fourier_df, folder, 'fourier_table', run_id, chunk_id, par_id, dry_run) logging.info(f"FFT data saved.") logging.info(f"FFT table contains {fourier_df.shape[0]} rows.") timing, full_df, _ = take_isrcs_with_inflection_last_week(fourier_df) save_fourier_s3_inflection_only(full_df, folder, 'series_with_inflection_last_week', run_id, chunk_id, par_id, dry_run) timing2, full_df_for_arima, _ = take_isrcs_with_inflection_last_week_GB(fourier_df) logging.info(f"{full_df_for_arima.shape[0]} rows loaded in {timing} seconds") # Arima timing, (cross_sec_df, ts_df), _ = arima_iscrs(full_df_for_arima) save_chunk_s3(cross_sec_df, folder, 'arima_cross', run_id, chunk_id, par_id, dry_run) save_chunk_s3(ts_df, folder, 'arima_timeseries', run_id, chunk_id, par_id, dry_run) return (cross_sec_df, ts_df) def get_argument(args, name, default_value=None): ''' Getting arguments from processing script ''' arg_name = "--" + name if arg_name in args: index = args.index("--" + name) + 1 return args[index] else: return default_value def get_mandatory_argument(args, name): ''' Getting mandatory arguments from processing script ''' res = get_argument(args, name) if res is None: raise "Missing script mandatory argument --{}".format(name) else: return res if __name__ == '__main__': # Dry run = True => data will NOT be saved DRY_RUN = get_argument(sys.argv, "dry-run", 'False') == 'True' # Run id RUN_ID = get_mandatory_argument(sys.argv, "run-id") # Folder FOLDER = get_mandatory_argument(sys.argv, "folder") # Splits the universe of ids into CHUNK_COUNT stable chunks CHUNK_COUNT = int(get_mandatory_argument(sys.argv, "chunk-count")) # We are running the script only for the below CHUNK_ID, out of CHUNK_COUNT chunks CHUNK_ID = int(get_mandatory_argument(sys.argv, "chunk-id")) # Data CSV File INPUT_FILE = 's3://dev-cucumbers/{}/MOMENTS_source_{}.csv'.format(FOLDER, RUN_ID) # Number of parallel processes to process the data PARALLELISM = 70 # int(get_argument(sys.argv, "parallelism", str(max(1, cpu_count() // 2)))) if DRY_RUN: logging.warning(f"DRY_RUN flag is True, no data will be saved !!!") logging.info(f"Processing Moments with run-id='{RUN_ID}'") logging.info(f"Reading data from CSV '{INPUT_FILE}'") logging.info(f"Computing chunk {CHUNK_ID} out of {CHUNK_COUNT} chunk(s)") logging.info(f"Data will be processed with parallelism={PARALLELISM}") # Reading input data file from S3 df = pd.read_csv(INPUT_FILE) # Fixing dates that are not read with the correct type in a CSV file df['ACTIVITY_DATE'] = pd.to_datetime(df['ACTIVITY_DATE']) df['ACTIVITY_DATE'] = df['ACTIVITY_DATE'].dt.date logging.info(f"Min activity date: {df['ACTIVITY_DATE'].min()}") logging.info(f"Max activity date: {df['ACTIVITY_DATE'].max()}") parallelism = PARALLELISM elapsed_time = split_work(df, parallelism, work_load, FOLDER, RUN_ID, CHUNK_ID, DRY_RUN) logging.info(f"Data processed in {elapsed_time} seconds")