# Dependencies # !pip -q install pyecharts absl-py # !pip -q install statsmodels # !pip -q install pmdarima import os import subprocess subprocess.check_call('pip install -r /opt/ml/processing/input/dependencies/requirements.txt', shell=True) import pandas as pd #import pyecharts as echarts import os 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.statespace.sarimax import SARIMAX from statsmodels.tsa.arima.model import ARIMA from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.seasonal import seasonal_decompose from statsmodels.tools.eval_measures import mse,rmse, meanabs from statsmodels.tsa.stattools import adfuller from statsmodels.tsa.statespace.tools import diff from scipy import fftpack from multiprocessing import Pool # Splits the universe of ids into CHUNK_COUNT stable chunks CHUNK_COUNT = 10 # We are running the script only for the below CHUNK_ID, out of CHUNK_COUNT chunks CHUNK_ID = 0 df = pd.read_csv('s3://dev-cucumbers/eimpara/Fourier/2023_data/data_for_2023_analysis_20231107-112155.csv') df['ACTIVITY_DATE'] = pd.to_datetime(df['ACTIVITY_DATE']) df['ACTIVITY_DATE'] = df['ACTIVITY_DATE'].dt.date print('Min activity date: ', df['ACTIVITY_DATE'].min()) print('Max activity date: ', df['ACTIVITY_DATE'].max()) def check_one_track_df(one_track_df): '''Checking that the length of the series has minimum 140 days''' return len(one_track_df)==140 def unique_isrc_df(isrc, df): '''Subsetting to have individual time series ahead of univariate time series analysis''' one_track_df = df[df['ISRC']==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', 'ACTIVITY_DATE', 'STREAMS']] return one_track_df.sort_values(by = 'ACTIVITY_DATE') def fourier_table(one_track_df, loop=60): '''Fast Fourier Transformation and returning a table.''' y = np.array(one_track_df['STREAMS']) y_fft_filtered = fftpack.fft(y).copy() freq = fftpack.fftfreq(len(y), d = 1/len(y)) cut_off = len(y)/loop y_fft_filtered[np.abs(freq) > cut_off]=0 Fourier = fftpack.ifft(y_fft_filtered) 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]] 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'] = np.array(one_track_df['ISRC']) return df_fourier def timing_val(func): def wrapper(*arg, **kw): '''From source: http://www.daniweb.com/code/snippet368.html''' t1 = time.time() res = func(*arg, **kw) t2 = time.time() return (t2 - t1), res, func.__name__ return wrapper @timing_val def compile_fourier_table(df): '''Assembling the new dataframe''' res = pd.DataFrame(columns=['ACTIVITY_DATE', 'STREAMS', 'Fourier', 'Inflection_Point', 'ISRC']) unique_isrcs = df['ISRC'].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: print('Error:',isrc,'data has length', len(one_track_df)) return res def save_dataframe_s3(df, chunk_id): '''Saving table to S3''' s3 = boto3.client('s3') bucket_name = 'dev-cucumbers' today = datetime.today().strftime('%Y%m%d-%H%M%S') filepath = "eimpara/Moments_2023_batches/Fourier_table_NEW_DATA_{}_{}.csv".format(chunk_id, today) csv_buffer = df.to_csv(index=False).encode('utf-8') s3.put_object(Body=csv_buffer, Bucket=bucket_name, Key=filepath) print(f"Table saved to S3 bucket: {bucket_name}, with file name: {filepath}") def data_for_chunk(df, chunk_id, chunk_count): IDs = df['ISRC'].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'].isin(ids)] return chunk_df chunk_table = data_for_chunk(df, CHUNK_ID, CHUNK_COUNT) fourier_table = compile_fourier_table(chunk_table) FourierTable = pd.DataFrame(fourier_table[1]) FourierTable['ACTIVITY_DATE'] = pd.to_datetime(FourierTable['ACTIVITY_DATE']).dt.date save_dataframe_s3(FourierTable, CHUNK_ID) def timing_val(func): def wrapper(*arg, **kw): t1 = time.time() res = func(*arg, **kw) t2 = time.time() return (t2 - t1), res, func.__name__ return wrapper @timing_val def load_full_data(): first_day_pred = df['ACTIVITY_DATE'].max() - timedelta(days=7) list_for_pred = FourierTable[(FourierTable['Inflection_Point']==1) & (FourierTable['ACTIVITY_DATE'] > first_day_pred)]['ISRC'].unique().tolist() subset_for_pred = df[df['ISRC'].isin(list_for_pred)].copy() return subset_for_pred def split_universe(df, n): def chunks(l, n): for i in range(0, len(l), n): yield l[i:i + n] isrcs = df['ISRC'].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): isrc_chunks = split_universe(full_df, n) splitted_df = [] for chunk in isrc_chunks: df_chunk = full_df[full_df['ISRC'].isin(chunk)] splitted_df.append(df_chunk) run_id = str(time.time() * 1000000) print("ISRCs have been splitted in {} dataframes for run id={}".format(n, run_id)) with Pool(n) as p: splitted_df_with_index = list(enumerate(splitted_df)) splitted_df_with_args = map(lambda x: (run_id, x[0], x[1]), splitted_df_with_index) p.map(f, splitted_df_with_args) print("Moments calculation finished.") # Save one chunk of work def save_chunk_s3(df, run_id, chunk_id): s3 = boto3.client('s3') bucket_name = 'dev-cucumbers' today = datetime.today().strftime('%Y%m%d-%H%M%S') filepath = "eimpara/Moments_2023_batches/ARIMA_chunks/{}_{}_{}.csv".format(run_id, chunk_id, today) csv_buffer = df.to_csv(index=False).encode('utf-8') s3.put_object(Body=csv_buffer, Bucket=bucket_name, Key=filepath) print(f"Table saved to S3 bucket: {bucket_name}, with file name: {filepath}") def auto_arima(df_isrc): '''Stepwise selection to decide autoregressive parameters''' 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', suppress_warnings=True, #stationary=False, stepwise=True) return auto_model def iterate_group_by_key(df, col_for_key, sort_data=True): if sort_data: df = df.sort_values(by=col_for_key) 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 arima_iscrs(df): '''ARIMA model, returns a table with model evaluation values''' new_df = pd.DataFrame(columns=['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']) count = 1 for subdf in iterate_group_by_key(df, ['ISRC']): df_isrc = subdf isrc = df_isrc['ISRC'].iloc[0] df_isrc = df_isrc.sort_values(by=['ACTIVITY_DATE']) try: cutoff_train_test_sets = 133 start_pred = 133 end_pred = 139 (train, test) = (df_isrc.iloc[:cutoff_train_test_sets], df_isrc.iloc[cutoff_train_test_sets:]) (end_test, end_train) = (len(test), len(train)) model = auto_arima(train) pred = model.predict(n_periods=7, return_conf_int=False) forecast_errors = np.subtract(np.array(test['STREAMS']), pred) 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' new_df = pd.concat([new_df, pd.DataFrame([{ 'ISRC': isrc, '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] }], columns=new_df.columns)]) if count % 100 ==0: print("Worker processed {} ISRCs.".format(count)) count = count+1 except Exception as e: print("Error while processing ISRC {}: '{}'".format(isrc, e)) return new_df def work_load(arg): (run_id, chunk_id, df) = arg print("Arguments {},{},{}".format(run_id, chunk_id, type(df))) timing, res, _ = arima_iscrs(df) save_chunk_s3(res, run_id, chunk_id) return res if __name__ == '__main__': parallelism = 10 timing, full_df, _ = load_full_data() print("{} rows loaded in {} seconds".format(full_df.shape[0], timing)) print("Splitting work into {} chunks".format(parallelism)) elapsed_time = split_work(full_df, parallelism, work_load) print("work done in {} seconds".format(elapsed_time))