<jupyter_start><jupyter_text>Network Traffic Forecasting (using time series data) In telco, accurate forecast of KPIs (e.g. network traffic, utilizations, user experience, etc.) for communication networks ( 2G/3G/4G/5G/wired) can help predict network failures, allocate resource, or save energy. In this notebook, we demonstrate a reference use case where we use the network traffic KPI(s) in the past to predict traffic KPI(s) in the future. We demostrate how to do **univariate forecasting** (predict only 1 series), and **multivariate forecasting** (predicts more than 1 series at the same time) using Project Chronos.For demonstration, we use the publicly available network traffic data repository maintained by the [WIDE project](http://mawi.wide.ad.jp/mawi/) and in particular, the network traffic traces aggregated every 2 hours (i.e. AverageRate in Mbps/Gbps and Total Bytes) in year 2018 and 2019 at the transit link of WIDE to the upstream ISP ([dataset link](http://mawi.wide.ad.jp/~agurim/dataset/)).   Helper functions This section defines some helper functions to be used in the following procedures. You can refer to it later when they're used.<jupyter_code>def get_drop_dates_and_len(df, target_col="total",allow_missing_num=3):
    """
    Find missing values and get records to drop
    """
    missing_num = df[target_col].isnull().astype(int).groupby(df[target_col].notnull().astype(int).cumsum()).sum()
    missing_num = missing_num[missing_num > 0]
    drop_datetimes = df.iloc[(missing_num.index + missing_num.cumsum() - missing_num).values].index
    drop_datetimes = drop_datetimes[missing_num > allow_missing_num]
    drop_missing_num = missing_num[missing_num > allow_missing_num]
    drop_len = drop_missing_num.values
    return drop_datetimes, drop_len
def rm_missing_weeks(start_dts, missing_lens, df):
    """
    Drop weeks that contains more than 3 consecutive missing values.
    If consecutive missing values across weeks, we remove all the weeks.
    """ 
    for start_time, l in zip(start_dts, missing_lens):
        start = start_time - pd.Timedelta(days=start_time.dayofweek)
        start = start.replace(hour=0, minute=0, second=0)
        start_week_end = start + pd.Timedelta(days=6)
        start_week_end = start_week_end.replace(hour=22, minute=0, second=0)

        end_time = start_time + l*pd.Timedelta(hours=2)
        if start_week_end < end_time:
            end = end_time + pd.Timedelta(days=6-end_time.dayofweek)
            end = end.replace(hour=22, minute=0, second=0)
        else:
            end = start_week_end
        df = df.drop(df[start:end].index)
    return df
def gen_dataset_matrix(dataset, look_back, target_col_indexes):
    """
    Generate input samples from rolling
    """
    X, Y = [], []
    if len(target_col_indexes) == 1:
        del_col_index = {0, 1}.difference(set(target_col_indexes)).pop()
        data = np.delete(dataset, del_col_index, axis=1)
    else:
        data = dataset
    for i in range(len(data) - look_back):
        X.append(data[i: (i + look_back)])
        Y.append(data[i + look_back, target_col_indexes])
    return np.array(X), np.array(Y)
def unscale(scaler, y, target_col_indexes):
    """
    data needs to be normalized (scaled) before feeding into models. 
    This is to inverse the effect of normlization to get reasonable forecast results.
    """
    dummy_feature_shape = scaler.scale_.shape[0]
    y_dummy = np.zeros((y.shape[0], dummy_feature_shape))
    y_dummy[:, target_col_indexes] = y
    y_unscale = scaler.inverse_transform(y_dummy)[:,target_col_indexes]
    return y_unscale
EPSILON = 1e-10
def sMAPE(y_true, y_pred, multioutput="uniform_average"):
    """
    Symmetric Mean Average Percentage Error
    """
    if multioutput not in ["raw_values", "uniform_average"]:
        raise ValueError("multioutput must be 'raw_values' or 'uniform_average', got {}".format(multioutput))
    output_errors = np.mean(100 * np.abs(y_true - y_pred)/(np.abs(y_true) + np.abs(y_pred) + EPSILON), axis=0,)
    if multioutput == "raw_values":
        return output_errors
    return np.mean(output_errors)
def plot_predict_actual_values(date, y_pred, y_test, ylabel):
    """
    plot the predicted values and actual values (for the test data)
    """
    fig, axs = plt.subplots(figsize=(16,6))

    axs.plot(date, y_pred,color='red', label='predicted values')
    axs.plot(date, y_test,color='blue', label='actual values')
    axs.set_title('the predicted values and actual values (for the test data)')

    plt.xlabel('test datetime')
    plt.ylabel(ylabel)
    plt.legend(loc='upper left')
    plt.show()<jupyter_output><empty_output><jupyter_text>Download raw dataset and load into dataframe Now we download the dataset and load it into a pandas dataframe. Steps are as below. * First, run the script `get_data.sh` to download the raw data. It will download the monthly aggregated traffic data in year 2018 and 2019 into `data` folder. The raw data contains aggregated network traffic (average MBPs and total bytes) as well as other metrics. * Second, run `extract_data.sh` to extract relavant traffic KPI's from raw data, i.e. `AvgRate` for average use rate, and `total` for total bytes. The script will extract the KPI's with timestamps into `data/data.csv`.* Finally, use pandas to load `data/data.csv` into a dataframe as shown below<jupyter_code>import os 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
raw_df = pd.read_csv("data/data.csv")<jupyter_output><empty_output><jupyter_text>Below are some example records of the data<jupyter_code>raw_df.head()<jupyter_output><empty_output><jupyter_text>Data pre-processing Now we need to do data cleaning and preprocessing on the raw data. Note that this part could vary for different dataset. For the network traffic data we're using, the processing contains 3 parts:1. Convert string datetime to TimeStamp2. Unify the measurement scale for `AvgRate` value - some uses Mbps, some uses Gbps 3. Handle missing data (fill or drop).<jupyter_code>df = pd.DataFrame(pd.to_datetime(raw_df.StartTime))
# we can find 'AvgRate' is of two scales: 'Mbps' and 'Gbps' 
raw_df.AvgRate.str[-4:].unique()
# Unify AvgRate value
df['AvgRate'] = raw_df.AvgRate.apply(lambda x: float(x[:-4]) if x.endswith("Mbps") else float(x[:-4]) * 1000)
df["total"] = raw_df["total"]
df.set_index("StartTime", inplace=True)
df.head()
full_idx = pd.date_range(start=df.index.min(), end=df.index.max(), freq='2H')
df = df.reindex(full_idx)
print("no. of n/a values:")
print(df.isna().sum())<jupyter_output>no. of n/a values:
AvgRate    3
total      3
dtype: int64<jupyter_text>Here, we drop weeks with more than 3 consecutive missing values and fill other missing values remained.<jupyter_code>drop_dts, drop_len = get_drop_dates_and_len(df)
df = rm_missing_weeks(drop_dts, drop_len, df)
df.ffill(inplace=True)<jupyter_output><empty_output><jupyter_text>Plot the data to see how the KPI's look like<jupyter_code>ax = df.plot(y='AvgRate',figsize=(16,6), title="AvgRate of network traffic data")
ax = df.plot(y='total',figsize=(16,6), title="total bytes of network traffic data")<jupyter_output><empty_output><jupyter_text>Feature Engineering For feature engineering, we use year, month, week, day of week and hour as features in addition to the target KPI values<jupyter_code>df["year"] = df.index.year
df["month"] = df.index.month
df["week"] = df.index.week
df["dayofweek"] = df.index.dayofweek
df["hour"] = df.index.hour
df.head()<jupyter_output><empty_output><jupyter_text>Data preparation Now we split the dataset into train and test.<jupyter_code>test_num = 1000
# we look back one week data which is of the frequency of 2h.
look_back = 84
test_split_index = test_num + look_back
target_dim = 1

#train_df
train_df = df[:-test_num]
test_df = df[-test_split_index:]
test_df = test_df.reset_index(drop=True)<jupyter_output><empty_output><jupyter_text>Then standardize train, test data and featues<jupyter_code>from sklearn.preprocessing import StandardScaler
standard_scaler = StandardScaler()
scaled_train = standard_scaler.fit_transform(train_df)
scaled_test = standard_scaler.transform(test_df)<jupyter_output><empty_output><jupyter_text>Last, we generate model input samples by sliding window along time axis. Since univariant and multivariant forcasting uses different input data shape, we prepare data for each type of forcasting respectively. * for _univariate_ forecasting, we forecast ``AvgRate`` only.<jupyter_code># for univeriate
# here we take "AvgRate" as an example
uni_target_col_indexes = [0]
uni_target_value = "AvgRate"

x_train, y_train = gen_dataset_matrix(scaled_train, look_back, uni_target_col_indexes)
x_test, y_test = gen_dataset_matrix(scaled_test, look_back, uni_target_col_indexes)
x_train.shape, y_train.shape, x_test.shape, y_test.shape<jupyter_output><empty_output><jupyter_text>* for _multivariate_ forecasting, we forecast `AvgRate` and `total` at the same time.<jupyter_code># for multivariate
multi_target_col_indexes = [0, 1]
multi_target_value = ["AvgRate","total"]

x_train_m, y_train_m = gen_dataset_matrix(scaled_train, look_back, multi_target_col_indexes)
x_test_m, y_test_m = gen_dataset_matrix(scaled_test, look_back, multi_target_col_indexes)
x_train_m.shape, y_train_m.shape, x_test_m.shape, y_test_m.shape<jupyter_output><empty_output><jupyter_text>Time series forecasting   Univariate forecasting For univariate forcasting, we use LSTMForecaster for forecasting.<jupyter_code>from zoo.chronos.model.forecast.lstm_forecaster import LSTMForecaster<jupyter_output>WARNING:tensorflow:From /home/junweid/intel_PR4/analytics-zoo/pyzoo/zoo/tfpark/zoo_optimizer.py:73: The name tf.train.Optimizer is deprecated. Please use tf.compat.v1.train.Optimizer instead.<jupyter_text>First we initiate a LSTMForecaster. * `feature_dim` should match the training data input feature, so we just use the last dimension of train data shape. * `target_dim` equals the variate num we want to predict. We set target_dim=1 for univariate forecasting.<jupyter_code># build model
lstm_config = {"lstm_units": [32] * 2, "lr":0.001}
forecaster = LSTMForecaster(target_dim=1, feature_dim=x_train.shape[-1], **lstm_config)<jupyter_output>WARNING:tensorflow:From /home/junweid/anaconda3/envs/automl/lib/python3.6/site-packages/tensorflow_core/python/ops/resource_variable_ops.py:1630: calling BaseResourceVariable.__init__ (from tensorflow.python.ops.resource_variable_ops) with constraint is deprecated and will be removed in a future version.
Instructions for updating:
If using Keras pass *_constraint arguments to layers.<jupyter_text>Then we use fit to train the model. Wait sometime for it to finish.<jupyter_code>%%time
forecaster.fit(x=x_train, y=y_train, batch_size=1024, epochs=50, distributed=False)<jupyter_output>WARNING:tensorflow:From /home/junweid/anaconda3/envs/automl/lib/python3.6/site-packages/tensorflow_core/python/ops/math_grad.py:1424: where (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.
Instructions for updating:
Use tf.where in 2.0, which has the same broadcast rule as np.where
Train on 7676 samples
Epoch 1/50
7676/7676 [==============================] - 4s 466us/sample - loss: 0.9872 - mean_squared_error: 0.9872
Epoch 2/50
7676/7676 [==============================] - 3s 359us/sample - loss: 0.9073 - mean_squared_error: 0.9073
Epoch 3/50
7676/7676 [==============================] - 3s 364us/sample - loss: 0.8573 - mean_squared_error: 0.8573
Epoch 4/50
7676/7676 [==============================] - 3s 356us/sample - loss: 0.8008 - mean_squared_error: 0.8008
Epoch 5/50
7676/7676 [==============================] - 3s 359us/sample - loss: 0.7405 - mean_squared_error: 0.7405
Epoch 6/50
7676/7676 [==============================] - 3s 362us/sampl[...]<jupyter_text>After training is finished. You can use the forecaster to do prediction and evaluation.<jupyter_code># make prediction
y_pred = forecaster.predict(x_test)<jupyter_output><empty_output><jupyter_text>Since we have used standard scaler to scale the input data (including the target values), we need to inverse the scaling on the predicted values too.<jupyter_code>y_pred_unscale = unscale(standard_scaler, y_pred, uni_target_col_indexes)
y_test_unscale = unscale(standard_scaler, y_test, uni_target_col_indexes)<jupyter_output><empty_output><jupyter_text>calculate the symetric mean absolute percentage error.<jupyter_code># evaluate with sMAPE
print("sMAPE is", sMAPE(y_test_unscale, y_pred_unscale))
# evaluate with mean_squared_error
from sklearn.metrics import mean_squared_error
print("mean_squared error is", mean_squared_error(y_test_unscale, y_pred_unscale))<jupyter_output>sMAPE is 9.510807442286207
mean_squared error is 6912.6167529738495<jupyter_text>multivariate forecasting  For multivariate forecasting, we use MTNetForecaster for forecasting.<jupyter_code>from zoo.chronos.model.forecast.mtnet_forecaster import MTNetForecaster<jupyter_output><empty_output><jupyter_text>First, we initialize a mtnet_forecaster according to input data shape. The lookback length is equal to `(long_series_num+1)*series_length` Details refer to [chronos docs](https://github.com/intel-analytics/analytics-zoo/tree/master/pyzoo/zoo/chronos).<jupyter_code>mtnet_forecaster = MTNetForecaster(target_dim=y_train_m.shape[-1],
                        feature_dim=x_train_m.shape[-1],
                        long_series_num=6,
                        series_length=12,
                        ar_window_size=6,
                        cnn_height=4
    )<jupyter_output>WARNING:tensorflow:From /home/junweid/anaconda3/envs/automl/lib/python3.6/site-packages/tensorflow_core/python/keras/initializers.py:94: calling TruncatedNormal.__init__ (from tensorflow.python.ops.init_ops) with dtype is deprecated and will be removed in a future version.
Instructions for updating:
Call initializer instance with the dtype argument instead of passing it to the constructor<jupyter_text>MTNet needs to preprocess the X into another format, so we call `MTNetForecaster.preprocess_input` on train_x and test_x.<jupyter_code># mtnet requires reshape of input x before feeding into model.
x_train_mtnet = mtnet_forecaster.preprocess_input(x_train_m)
x_test_mtnet = mtnet_forecaster.preprocess_input(x_test_m)<jupyter_output><empty_output><jupyter_text>Now we train the model and wait till it finished.<jupyter_code>%%time
hist = mtnet_forecaster.fit(x = x_train_mtnet, y = y_train, batch_size=1024, epochs=20)<jupyter_output>Train on 7676 samples
Epoch 1/20
7676/7676 [==============================] - 18s 2ms/sample - loss: 0.9270 - mean_squared_error: 1.3800
Epoch 2/20
7676/7676 [==============================] - 3s 447us/sample - loss: 0.8441 - mean_squared_error: 1.2350
Epoch 3/20
7676/7676 [==============================] - 3s 447us/sample - loss: 0.7853 - mean_squared_error: 1.1383
Epoch 4/20
7676/7676 [==============================] - 3s 446us/sample - loss: 0.7286 - mean_squared_error: 0.9931
Epoch 5/20
7676/7676 [==============================] - 3s 447us/sample - loss: 0.6647 - mean_squared_error: 0.8576
Epoch 6/20
7676/7676 [==============================] - 3s 446us/sample - loss: 0.5952 - mean_squared_error: 0.7253
Epoch 7/20
7676/7676 [==============================] - 3s 446us/sample - loss: 0.5420 - mean_squared_error: 0.5940
Epoch 8/20
7676/7676 [==============================] - 3s 447us/sample - loss: 0.5082 - mean_squared_error: 0.5153
Epoch 9/20
7676/7676 [=============================[...]<jupyter_text>Use the model for prediction and inverse the scaling of the prediction results<jupyter_code>y_pred_m = mtnet_forecaster.predict(x_test_mtnet)
y_pred_m_unscale = unscale(standard_scaler, y_pred_m, multi_target_col_indexes)
y_test_m_unscale = unscale(standard_scaler, y_test_m, multi_target_col_indexes)
# evaluate with sMAPE
print("sMAPE is",sMAPE(y_pred_m_unscale, y_test_m_unscale, multioutput="raw_values"))
# evaluate with mean_squared_error
print("mean squared error is", mean_squared_error(y_test_m_unscale, y_pred_m_unscale, multioutput="raw_values"))<jupyter_output>sMAPE is [10.22883616 11.34957483]
mean squared error is [9.44256767e+03 1.04967574e+22]<jupyter_text>plot actual and prediction values for `AvgRate` KPI<jupyter_code>test_date=df[-test_num:].index
plot_predict_actual_values(test_date, y_pred_m_unscale[:,0], y_test_m_unscale[:,0], ylabel=multi_target_value[0])<jupyter_output><empty_output><jupyter_text>plot actual and prediction values for `total bytes` KPI<jupyter_code>plot_predict_actual_values(test_date, y_pred_m_unscale[:,1], y_test_m_unscale[:,1], ylabel=multi_target_value[1])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12410158121705797, 'fraction_numerical': 0.03497843794920939, 'mean_word_length': 5.522265625, 'pattern_counts': {'":': 3, '<': 90, '<?xml version=': 0, '>': 92, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '118050', 'license': 'Apache-2.0', 'n_tokens_mistral': 6041, 'n_tokens_neox': 5486, 'n_words': 1438, 'path': 'pyzoo/zoo/chronos/use-case/network_traffic/network_traffic_model_forecasting.ipynb', 'repo_name': 'block4chain/analytics-zoo', 'size': '16696'}