<jupyter_start><jupyter_text>Understanding Anomalies in Time-Series - Is something going on?<jupyter_code>!pip install burst_detection
!pip install torch
!pip install torchvision
!apt-get install subversion
!svn export https://github.com/FreelancerFel/Anomalies_Detection/trunk/datasets

# Source for burst_detection implementation
# https://github.com/nmarinsek/burst_detection
# Source for Numenta Anomaly Benchmark
# https://github.com/numenta/NAB<jupyter_output>Requirement already satisfied: burst_detection in /usr/local/lib/python3.6/dist-packages (0.1.3)
Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from burst_detection) (1.18.5)
Requirement already satisfied: pandas in /usr/local/lib/python3.6/dist-packages (from burst_detection) (1.1.4)
Requirement already satisfied: mpmath in /usr/local/lib/python3.6/dist-packages (from burst_detection) (1.1.0)
Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas->burst_detection) (2018.9)
Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas->burst_detection) (2.8.1)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.7.3->pandas->burst_detection) (1.15.0)
Requirement already satisfied: torch in /usr/local/lib/python3.6/dist-packages (1.7.0+cu101)
Requirement already satisfied: typing-extensions in /usr/local/lib/pyt[...]<jupyter_text>Objective* Understand anomalies in time-series* Ability to implement anomalies detection using forecasting with ARIMA, Kleinburg's Algorithm for Burst Detection, and Pytorch implementation of Autoencoder.<jupyter_code>import warnings
warnings.filterwarnings('ignore')
# load packages

import pandas as pd
import numpy as np
%pylab inline 
### Dataset to play with?
from itertools import accumulate
from plotnine import *
import burst_detection as bd

from sklearn import preprocessing

import math
import statsmodels.api as sm
import statsmodels.tsa.api as smt
from sklearn.metrics import mean_squared_error
from matplotlib import pyplot
import matplotlib.pyplot as plt<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>Data We use twitter data from the Numenta Anomaly Benchmark (NAB) v1.1. NAB is a novel benchmark for evaluating algorithms for anomaly detection in streaming, real-time applications. The twitter data contain number of tweets related to ten company ranging from Apple, Amazon to Google between the period of Feb 2015 to April 2015<jupyter_code>### Load All twitter data, identify where is the 'unusual traffic'
tweets = ['AAPL','AMZN','CRM','CVS','FB','GOOG','IBM','KO','PFE','UPS']
df = pd.read_csv('./datasets/Twitter_volume_AAPL.csv', parse_dates=True)
df['key'] = 'AAPL'
for i in tweets[1:]:
    filename = './datasets/Twitter_volume_'+str(i)+'.csv'
    file = pd.read_csv(filename, parse_dates=True)
    file['key'] = i
    df = pd.concat([df,file], sort=False)
df<jupyter_output><empty_output><jupyter_text>The dataset is collected every five minutes, but we are not interested in that small of an inteval for time period so it is sum at hour level.Additionally, we also create a column for aggregate tweets for that company up until that moment and total number of tweets (for all 10 companies) for that hours. The reasoning for the total is to create a 'normal' behaviour for benchmarking<jupyter_code>df['timestamp'] =  pd.to_datetime(df['timestamp'], format='%Y-%m-%d %H:%M:%S')
df['timestamp'] = df['timestamp'].apply(lambda dt: datetime.datetime(dt.year, dt.month, dt.day, dt.hour,dt.minute//60))
df_alt = df.groupby(['timestamp'])['value'].agg('sum').reset_index(inplace = False)
df = df.groupby(['timestamp','key'])['value'].agg('sum').reset_index(inplace = False)
df = df.sort_values(by=['key','timestamp'])
original = []
for i in tweets:
    output = [*accumulate(list(df[df['key']==i]['value']), lambda a, b: a+b)]
    original= original + output
df['aggregate'] = original
df = df[df['value']!=0]
df = df.merge(df_alt, on=['timestamp'], how='left')
df.columns = ['timestamp','key','value','aggregate','total']
df.head(n=20)<jupyter_output><empty_output><jupyter_text>In NAB, we also have the answer. lA json file labelling all 'anomalies' datapoint which is given to us to help with validation.<jupyter_code>### Read Answer Data
import json
# read file
with open('./datasets/combined_labels.json', 'r') as myfile:
    data=myfile.read()
# parse file
obj = json.loads(data)
datetimenew = []
keyword = []
for i in tweets:
    filename = 'realTweets/Twitter_volume_'+str(i)+'.csv'
    size = len(obj[filename])
    datetimenew += obj[filename]
    keyword += [i]*size
target = pd.DataFrame({'key':keyword,'timestamp':datetimenew,'outlier':1})
target['timestamp'] =  pd.to_datetime(target['timestamp'], format='%Y-%m-%d %H:%M:%S')
target['timestamp'] = target['timestamp'].apply(lambda dt: datetime.datetime(dt.year, dt.month, dt.day, dt.hour,dt.minute//60))
target.head()
df = df.merge(target, on=['key','timestamp'], how='left')
df = df.fillna(0)<jupyter_output><empty_output><jupyter_text>We now have a dataset with six columns containing everything we need to run some basic anomaly/outlier detection.<jupyter_code>df.head()<jupyter_output><empty_output><jupyter_text>Visualization Before we get into the detail, let's see what is anomalies in NAB dataset and how does it look like. Some basic plot of each time-series and anomalies are shown below.<jupyter_code>ggplot(df, aes(x='timestamp',y='value', group='key', color='key'))+ geom_line()+ theme(
  axis_text_x = element_text(size = 5.5,angle = 90, hjust = 20, vjust = -400)
)+geom_point(df[df['outlier']==1],aes(x='timestamp',y='value', group='key'), color='black')+facet_wrap('~key', scales='free')<jupyter_output><empty_output><jupyter_text>What you may also noticed is that anomalies point that were given by NAB benchmark wasn't entirely accurate in all case or at the very least up for interpretation.<jupyter_code>obj['realTweets/Twitter_volume_AAPL.csv']
ggplot(df[df['key']=='AAPL'], aes(x='timestamp',y='value', group='key', color='key'))+ geom_line()+ theme(
  axis_text_x = element_text(size = 5.5,angle = 90, hjust = 20, vjust = -400)
)+geom_point(df[(df['outlier']==1)&(df['key']=='AAPL')],aes(x='timestamp',y='value', group='key'), color='black')<jupyter_output><empty_output><jupyter_text>We also run visualization by 'aggregate data', this is less clear.<jupyter_code>ggplot(df, aes(x='timestamp',y='aggregate', group='key', color='key'))+ geom_line()+ theme(
  axis_text_x = element_text(size = 5.5,angle = 90, hjust = 20, vjust = -400)
)+geom_point(df[df['outlier']==1],aes(x='timestamp',y='aggregate', group='key'), color='black')+facet_wrap('~key', scales='free') 
ggplot(df[df['key']=='AAPL'], aes(x='timestamp',y='aggregate', group='key', color='key'))+ geom_line()+ theme(
  axis_text_x = element_text(size = 5.5,angle = 90, hjust = 20, vjust = -400)
)+geom_point(df[(df['outlier']==1)&(df['key']=='AAPL')],aes(x='timestamp',y='aggregate', group='key'), color='black')<jupyter_output><empty_output><jupyter_text>Anomalies After some basic visualization to understand what kind of anomalies existed in this example dataset., we now just need to give a brief explanation of an anomaly (or outlier). It is the deviation in a quantity from its expected value. Anomaly detection is the search of that which ‘does not conform’ with expected value. * Point Anomaly![Point Anomaly](https://github.com/FreelancerFel/Anomalies_Detection/blob/master/global-outlier.png?raw=1)* Contextual Anomaly![Anomaly](https://github.com/FreelancerFel/Anomalies_Detection/blob/master/contextual-outlier.png?raw=1)* Collective Anomaly![Anomaly](https://github.com/FreelancerFel/Anomalies_Detection/blob/master/collective-outliers.png?raw=1)  Looking into the future - Using Forecasting to detect anomalies One of the most simple way of detecting anomalies is to use time-series prediction. Essentially, anomalies is about finding datapoint that is out of ordinary. Time series forecasting helps us estimate that ordinary which we can then use to detect anomalies by comparing them with actual data.**It does requires you to set a threshold, but this can be navigate with multiple threshold if you are unsure.**This obviously relies on good forecasting model, but it also means that any forecasting model is compatible with this idea.* Exponential Smoothing* ARIMA* XGBoost - Lag FeatureIn this instance, we will simply use ARIMA which is essentially lag feature regression.<jupyter_code>AnomaliesData = df[df['key']=='AAPL']
train, test = (AnomaliesData.value[0:-48]), (AnomaliesData.value[-48:])
train_log, test_log = np.log10(train), np.log10(test) # take log to scale
AnomaliesData.head()<jupyter_output><empty_output><jupyter_text>Determine P,D,Q through ACF and PACF table, alternatively use grid search according to either AIC, AICc or BIC value<jupyter_code>fig, axes = plt.subplots(1, 2, figsize=(15,4))
fig = sm.graphics.tsa.plot_acf(train, lags=40, ax=axes[0])
fig = sm.graphics.tsa.plot_pacf(train, lags=40, ax=axes[1])
order_pdq = (1, 1, 1)
seasonal_order_pdq = (0, 1, 1, 6)
datapoint = [x for x in train_log]
model = sm.tsa.SARIMAX(datapoint, order=order_pdq, seasonal_order=seasonal_order_pdq,enforce_stationarity=False,enforce_invertibility=False)
res = model.fit(disp=False)
print(res.summary())<jupyter_output>Statespace Model Results                                
=========================================================================================
Dep. Variable:                                 y   No. Observations:                 1276
Model:             SARIMAX(1, 1, 1)x(0, 1, 1, 6)   Log Likelihood                 181.208
Date:                           Sun, 08 Nov 2020   AIC                           -354.416
Time:                                   23:36:19   BIC                           -333.857
Sample:                                        0   HQIC                          -346.691
                                          - 1276                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------[...]<jupyter_text>We are not trying to predict everything from scratch, but rather use information at each stage to predict what is likely to be the next stage. So instead of train/test split that are norm in forecasting, we reuse the data at every stage and make prediction only for the next datapoint.<jupyter_code>predictions = list()
predict_log = list()
for t in range(len(test_log)):
    model = sm.tsa.SARIMAX(datapoint, order=order_pdq, seasonal_order=seasonal_order_pdq,enforce_stationarity=False,enforce_invertibility=False)
    model_fit = model.fit(disp=0)
    output = model_fit.forecast()
    predict_log.append(output[0])
    yhat = 10**output[0]
    predictions.append(yhat)
    obs = test_log.iloc[t]
    datapoint.append(obs)
predictions = pd.Series(predictions)
predictions.index = test.index
figsize=(15, 7)
plt.figure(figsize=figsize)
pyplot.plot(train,label='train', color='blue')
pyplot.plot(test,label='test', color='purple')
pyplot.plot(predictions, color='green',label='Predicted')
pyplot.legend(loc='upper right')
pyplot.show()<jupyter_output><empty_output><jupyter_text>Obviously, if you want to look at anomalies in a certain area, you shouldn't train your forecast model using that data. So we reserved some data at the end for some test.<jupyter_code>figsize=(15, 7)
plt.figure(figsize=figsize)
pyplot.plot(test,label='test', color='purple')
pyplot.plot(predictions, color='green',label='Predicted')
pyplot.legend(loc='upper right')
pyplot.show()
predicted_df=pd.DataFrame()
predicted_df['load_date']=AnomaliesData['timestamp'][-48:]
predicted_df['actuals']=test
predicted_df['predicted']=predictions
predicted_df.reset_index(inplace=True)
del predicted_df['index']
predicted_df.head()<jupyter_output><empty_output><jupyter_text>What we have now is the real and predicted data, we just need to look at the average and distribution of the 'error' of our model to create potential approximate 'normal' pattern.Subsequently, anything outside of that pattern is a potential anomalies.<jupyter_code>import numpy as np
def anomalies(df,rollingwindows):
    df['error']=df['actuals']-df['predicted']
    df['meanerror'] = df['error'].rolling(window=rollingwindows).mean()
    df['deviationerror'] = df['error'].rolling(window=rollingwindows).std()
    df['-3sderror'] = df['meanerror'] - (3 * df['deviationerror'])
    df['3sderror'] = df['meanerror'] + (3 * df['deviationerror'])
    df['-2sderror'] = df['meanerror'] - (2 * df['deviationerror'])
    df['2sderror'] = df['meanerror'] + (2 * df['deviationerror'])
    df['-1sderror'] = df['meanerror'] - (1 * df['deviationerror'])
    df['1sderror'] = df['meanerror'] + (1 * df['deviationerror'])
    cut_list = df[['error', '-3sderror', '-2sderror', '-1sderror', 'meanerror', '1sderror', '2sderror', '3sderror']]
    cut_sort = np.sort(cut_list.values)
    df['Anomalies_Strength'] = [(lambda number: np.where(cut_sort == df['error'][number])[1][0])(number)\
                                for number in range(len(df['error']))]
    severity = {0: 'Severe Anomaly', 1: 'Anomaly', 2: 'Unusual', 3: 'Normal', 4: 'Normal',
                5: 'Unusual', 6: 'Anomaly', 7: 'Severe Anomaly'}
    df['Anomalies_Strength'] =  df['Anomalies_Strength'].map(severity)
    df['anomaly_points'] = np.where(df['Anomalies_Strength'] == 'Severe Anomaly', df['error'], np.nan)
    df = df.sort_values(by='load_date', ascending=False)
    df.load_date = pd.to_datetime(df['load_date'].astype(str), format="%Y-%m-%d")
    df[-rollingwindows:]['anomaly_points'] = np.NaN
    df[-rollingwindows:]['Anomalies_Strength'] = 'Normal'
    df['upper_boundary'] = df['predicted']+df['1sderror']
    df['lower_boundary'] = df['predicted']+df['-1sderror']
    df['upper_boundary_high'] = df['predicted']+df['3sderror']
    df['lower_boundary_high'] = df['predicted']+df['-3sderror']
    df = df.drop(columns=['1sderror', '-1sderror','2sderror','-2sderror','-3sderror','3sderror'])
    return df
classify_df=anomalies(predicted_df,12)
ggplot(classify_df)+ geom_line(aes(x='load_date',y='actuals')) \
+ geom_point(aes(x='load_date',y='actuals',color='Anomalies_Strength'))+\
  scale_color_manual(values = {'Normal':"black",'Anomaly':"orange",'Severe Anomaly':"red",'Unusual':'yellow'})\
+theme(axis_text_x = element_text(size = 5.5,angle = 90, hjust = 20, vjust = -400))\
+geom_line(aes(x='load_date',y='upper_boundary'),color='orange',linetype="dotted")\
+geom_line(aes(x='load_date',y='lower_boundary'),color='orange',linetype="dotted")\
+geom_line(aes(x='load_date',y='upper_boundary_high'),color='red',linetype="dotted")\
+geom_line(aes(x='load_date',y='lower_boundary_high'),color='red',linetype="dotted")
classify_df.head()<jupyter_output><empty_output><jupyter_text>Looking into the past - Kleinberg's Burst Detection http://www.cs.cornell.edu/home/kleinber/bhs.pdfKleinberg’s algorithm identifies periods of time in which an event of particular interest is uncharacteristically frequent, or “bursty.”  It can be used to detect bursts in a continuous stream or batches of events by assumeing that there are only two possible states, the baseline state and the 'bursty' state. The probability of event of interest happening in baseline state is equal to the total number of event of interest divided by total number of events. The 'bursty' state probability is constant multiplied by baseline probability. * pybursts (Continuous stream) = https://github.com/romain-fontugne/pybursts (Forked version which fix some problem)* burst_detection (Batch event) = https://github.com/nmarinsek/burst_detection The total cost of transitioning from normal state to burst is equal to the sum of the two functions below. With the cost found, we can then find the optimal state sequence that tell us what is the state at each timepoint by measuring if the new state has been reached or not. The goodness of fit between the observed proportion and the expected probability of each state.$\sigma= - ln[\binom{d_t}{r_t}p_i^{rt}(1-p_i)^{d_t-r_t}]$<jupyter_code>import sympy

def fitcost(normal,interest,p):
    return -np.log(np.float(sympy.binomial(normal,interest))*(p**interest)*((1-p)**(normal-interest)))<jupyter_output><empty_output><jupyter_text>The difficulty of transitioning from the previous state to the next state.$\tau= - (i_{next}-i_{prev})*\gamma*ln(n)$<jupyter_code>def transitioncost(previousstate,nextstate,gamma,n):
    if previousstate>=nextstate:
        return 0
    else: 
        return (nextstate-previousstate) * gamma * np.log(n)
def burstdetection(interest,normal,totaltimepoint,activestatecost,gamma,smooth):
    
    k = 2 
    # smooth time_series, can be ignored if not needed.
    if smooth > 1:
        temp_p = interest/normal
        temp_p = temp_p.rolling(window=smooth, center=True).mean()
        interest = temp_p*normal
        real_n = sum(~np.isnan(interest)) 
    else: 
        real_n = totaltimepoint
          
    #calculate the expected proportions for states 0 and 1
    p = {}
    p[0] = np.nansum(interest) / float(np.nansum(normal)) #baseline state
    p[1] = p[0] * activestatecost  #cost to move from normal state to a burst state. 
    if p[1] > 1: #active state can't be more than maximum of proportuib
      p[1] = 0.99999

    #initialize matrices to hold the costs and optimal state sequence
    cost = np.zeros([totaltimepoint,k])
    q = np.zeros([totaltimepoint,1])

    for t in range(int((smooth-1)/2),(int((smooth-1)/2))+real_n):

        for j in range(k): # There are two states, which one has less cost: 'Active' or Normal
            if t == int((smooth-1)/2):
                cost[t,j] = fitcost(normal[t],interest[t],p[j])
            else:
                cost[t,j] = transitioncost(q[t-1],j,gamma,real_n) + fitcost(normal[t],interest[t],p[j])

        #add the state with the minimum cost to the optimal state sequence
        q[t] = np.where(cost[t,:] == min(cost[t,:]))

    return q,p
AnomaliesData_Start = df[df['key']=='AAPL']
min_max_scaler = preprocessing.MinMaxScaler()
np_scaled = min_max_scaler.fit_transform(AnomaliesData_Start[['value','total']])
AnomaliesData = pd.DataFrame(np_scaled*100, columns = ['value','total']) 
#Scaled the data because of issue with binomial coefficient
r = pd.Series(AnomaliesData['value']).astype('float')
d = pd.Series(AnomaliesData['total']).astype('float')<jupyter_output><empty_output><jupyter_text>Burst Detection Function: burst_detection(r,d,len(r),s,d,smooth_win)* r = number of events of interested* d = number of all events* s = distance between normal state and a burst state* gamma = cost to move from normal state to a burst state.* smooth = smoothing windows<jupyter_code>[q,p] = burstdetection(r,d,len(r),2,1,smooth=3)
bursts = bd.enumerate_bursts(q, 's='+str(3)+', g='+str(1))
bursts = bd.burst_weights(bursts, r, d, p)
data_r = pd.DataFrame({'day':range(len(r)),'data':r,'total':d,'outlier':AnomaliesData_Start['outlier']})
data_r['Prediction'] = 0
for index, row in bursts.iterrows():
    data_r['Prediction'] = np.where(data_r['day']>=row.begin,np.where(data_r['day']<=row.end,1,data_r['Prediction']), data_r['Prediction'])
ggplot(data_r, aes(x='day',y='data'))+ geom_line()+ theme(
  axis_text_x = element_text(size = 5.5,angle = 90, hjust = 20, vjust = -400)
)+geom_point(data_r[(data_r['outlier']==1)],aes(x='day',y='data'), color='red',shape = 7)+geom_point(data_r[(data_r['Prediction']==1)],aes(x='day',y='data'), color='blue',shape = 6)<jupyter_output><empty_output><jupyter_text>The strength, or weight, of a burst (that begins at time point t1 and ends at time point t2) can be estimated with the following function:<jupyter_code>def plot_burst_timeline(bursts, timepoints):
    
    f, ax = plt.subplots(figsize=(8, 0.5))
    ax.set(xlim=(0, timepoints), ylabel="", xlabel="")

    #create boxes around bursting periods
    for index, burst in bursts.iterrows():

        #define outline positions
        y = 0.25
        xstart = burst['begin'] -1 
        width = burst['end'] - burst['begin'] +1

        ax.add_patch(Rectangle((xstart,y), width, height=0.5, 
                               facecolor='#00bbcc', edgecolor='none', linewidth=1))

    ax.xaxis.set_visible(False)

    plt.axhline(0.5, linewidth=1, color='k', alpha=1, zorder=0.5)
    plt.show()
plot_burst_timeline(bursts, len(r))<jupyter_output><empty_output><jupyter_text>Finding the Pattern - CNN Autoencoders  Last, but by no mean least. There might be some situation where you lack understanding of the data to be able to correctly estimated boundary or level of energy for either forecasting or burst detection.Essentially, Autoencoder try to reconstruct the data based on what is learn before. If the autoencoder is good enough, it should be able to generalized to most situation and retain most of the standard information of time-series.<jupyter_code>import time
import numpy as np
import torch
import torch.nn.functional as F
from torchvision import datasets
from torchvision import transforms
from torch.utils.data import DataLoader
from sklearn.model_selection import train_test_split

if torch.cuda.is_available():
    torch.backends.cudnn.deterministic = True
    
def get_value_from_df(df):
    return df.value.to_list()

def normalize(values,mean=None,std=None):
    if mean is None:
        mean = np.mean(values)
    values -= mean
    if std is None:
        std = np.std(values)
    values /= std
    return values, mean, std

def create_sequences(values, time_steps=24):
    output = []
    for i in range(len(values) - time_steps):
        output.append(values[i : (i + time_steps)])
    # Convert 2D sequences into 3D as we will be feeding this into
    # a convolutional layer.
    return np.expand_dims(output, axis=2)<jupyter_output><empty_output><jupyter_text>Using the same data and splitting into test and train, we then decided to use 24 hours of data for each autoencoders. Hence, 24 timesteps.<jupyter_code>Autoencoder_Data = df[df['key']=='AAPL']
# Get the `value` column from the training dataframe.
training_value = get_value_from_df(Autoencoder_Data)
# Normalize `value` and save the mean and std we get
training_value, training_mean, training_std = normalize(training_value)
TIME_STEPS = 24
batch_size = 128
train = create_sequences(training_value,24)
train, test = (train[0:1000,:,:], train[1000:,:,:])
X_train = torch.FloatTensor([train]).cuda().reshape(-1, 24, 1)
X_test = torch.FloatTensor([test]).cuda().reshape(-1, 24, 1)
print("training input shape: ", X_train.shape)
print("testing input shape: ", X_test.shape)
train_loader = DataLoader(dataset=X_train, 
                          batch_size=batch_size, 
                          shuffle=True)

test_loader = DataLoader(dataset=X_test, 
                         batch_size=batch_size, 
                         shuffle=False)
##########################
### SETTINGS
##########################

# Device
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
print('Device:', device)

# Hyperparameters
learning_rate = 0.001
num_epochs = 100<jupyter_output>Device: cuda:0<jupyter_text>Autoencoder model is created here. The model will take input and return output of the same shape.<jupyter_code>class ConvolutionalAutoencoder(torch.nn.Module):

    def __init__(self):
        super(ConvolutionalAutoencoder, self).__init__()
        

        ### ENCODER
        self.conv_1 = torch.nn.Conv1d(in_channels=1,
                                      out_channels=32,
                                      kernel_size=7,
                                      stride=2,
                                      padding=1)        
        self.dropout_1 = torch.nn.Dropout(p=0.2)                                 
        self.conv_2 = torch.nn.Conv1d(in_channels=32,
                                      out_channels=16,
                                      kernel_size=7,
                                      stride=2,
                                      padding=1)                 
        self.dropout_2 = torch.nn.Dropout(p=0.2)       

        ### DECODER                       
        self.deconv_1 = torch.nn.ConvTranspose1d(in_channels=16,
                                                 out_channels=32,
                                                 kernel_size=7,
                                                 stride=2,
                                                 padding=1)
        self.dropout_3 = torch.nn.Dropout(p=0.2)  
        self.deconv_2 = torch.nn.ConvTranspose1d(in_channels=32,
                                                 out_channels=32,
                                                 kernel_size=7,
                                                 stride=2,
                                                 padding=1)
        self.deconv_3 = torch.nn.ConvTranspose1d(in_channels=32,
                                                 out_channels=1,
                                                 kernel_size=7,
                                                 stride=1,
                                                 padding=1)
        
        
    def forward(self, x):
        ### Input: (batch_size, seq_len, input_size)
        ### ENCODER
        batch_size = x.size()[0]
        seq_length = x.size()[1]
        x = x.view(batch_size,1,seq_length)
        x = self.conv_1(x)
        x = F.relu(x)
        x = self.dropout_1(x)
        x = self.conv_2(x)
        x = F.relu(x)
        x = self.dropout_2(x)
        
        ### DECODER
        x = self.deconv_1(x)
        x = F.relu(x)
        x = self.dropout_3(x)
        x = self.deconv_2(x)
        x = F.relu(x)
        x = self.deconv_3(x)
        x = x[:,:,:seq_length]
        x = x.view(batch_size,seq_length,1)
        return x

    
# torch.manual_seed(random_seed)
model = ConvolutionalAutoencoder()
model = model.to(device)

optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)  
model<jupyter_output><empty_output><jupyter_text>Train the model.<jupyter_code>start_time = time.time()
for epoch in range(num_epochs):
    for batch_idx, (features) in enumerate(train_loader):
        
        # don't need labels, only the images (features)
        features = features.to(device)

        ### FORWARD AND BACK PROP
        logits = model(features)
        cost = F.mse_loss(logits, features)
        optimizer.zero_grad()
        
        cost.backward()
        
        ### UPDATE MODEL PARAMETERS
        optimizer.step()
        
        ### LOGGING
        if not batch_idx % 5 and not epoch % 5:
            print ('Epoch: %03d/%03d | Batch %03d/%03d | Cost: %.4f' 
                   %(epoch+1, num_epochs, batch_idx, 
                     len(X_train)//batch_size, cost))
            
    print('Time elapsed: %.2f min' % ((time.time() - start_time)/60))
    
# print('Total Training Time: %.2f min' % ((time.time() - start_time)/60))<jupyter_output>Epoch: 001/100 | Batch 000/007 | Cost: 0.7682
Epoch: 001/100 | Batch 005/007 | Cost: 0.8018
Time elapsed: 0.00 min
Time elapsed: 0.00 min
Time elapsed: 0.00 min
Time elapsed: 0.00 min
Time elapsed: 0.00 min
Epoch: 006/100 | Batch 000/007 | Cost: 0.3898
Epoch: 006/100 | Batch 005/007 | Cost: 0.2225
Time elapsed: 0.00 min
Time elapsed: 0.00 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Epoch: 011/100 | Batch 000/007 | Cost: 0.3602
Epoch: 011/100 | Batch 005/007 | Cost: 0.1763
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Epoch: 016/100 | Batch 000/007 | Cost: 0.3334
Epoch: 016/100 | Batch 005/007 | Cost: 0.3780
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Epoch: 021/100 | Batch 000/007 | Cost: 0.2987
Epoch: 021/100 | Batch 005/007 | Cost: 0.1593
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapsed: 0.01 min
Time elapse[...]<jupyter_text>Detecting anomaliesAnomalies is detected by comparing input and output and determine those with mean absolute error that is higher than threshold.If the MEA loss is greater than  threshold value then it can be infer that the model is seeing a pattern that it isn't familiar with which is a potential anomaly. In this case, we used maximum error from training set.<jupyter_code>x_train_pred = model(X_train)
x_train_pred = x_train_pred.detach().cpu().numpy() 
train_mae_loss = np.mean(np.abs(x_train_pred - train), axis=1)

plt.hist(train_mae_loss, bins=50)
plt.xlabel("Train MAE loss")
plt.ylabel("No of samples")
plt.show()

# Get reconstruction loss threshold.
threshold = np.max(train_mae_loss)
print("Reconstruction error threshold: ", threshold)<jupyter_output><empty_output><jupyter_text>We now validate it by looking at test set for index of time that has high MAE.<jupyter_code>x_test_pred = model(X_test)
x_test_pred = x_test_pred.detach().cpu().numpy() 
test_mae_loss = np.mean(np.abs(x_test_pred - test), axis=1)
test_mae_loss = test_mae_loss.reshape((-1))

plt.hist(test_mae_loss, bins=50)
plt.xlabel("test MAE loss")
plt.ylabel("No of samples")
plt.show()

# Detect all the samples which are anomalies.
anomalies = (test_mae_loss > threshold).tolist()
index = np.where(anomalies)
print("Number of anomaly samples: ", np.sum(anomalies))
print("Indices of anomaly samples: ", np.where(anomalies))
Test_Autoencoder = Autoencoder_Data[1000:]
anomalous_data_indices = []
for data_idx in range(24 - 1, len(Test_Autoencoder) - 24 + 1):
    time_series = range(data_idx - 24 + 1, data_idx)
    if any([anomalies[j] for j in time_series]):
        Test_Autoencoder['outlier'].iloc[data_idx,:] = 1
        anomalous_data_indices.append(data_idx)
Test_Autoencoder[Test_Autoencoder['outlier']==1]
df_subset = Test_Autoencoder.iloc[anomalous_data_indices, :]
plt.subplots_adjust(bottom=0.2)
plt.xticks(rotation=25)
ax = plt.gca()

values = pd.Series(Test_Autoencoder["value"])
plt.plot(values, label="test data")
values = pd.Series(df_subset["value"])
plt.plot(values, label="anomalies", color="r")

plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.137189811010682, 'fraction_numerical': 0.024157764995891538, 'mean_word_length': 3.7960277427490543, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 106, 'https://': 8, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 44, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '471768', 'license': 'Apache-2.0', 'n_tokens_mistral': 10423, 'n_tokens_neox': 9636, 'n_words': 2481, 'path': 'Anomalies_Detection.ipynb', 'repo_name': 'FreelancerFel/Anomalies_Detection', 'size': '30425'}