<jupyter_start><jupyter_text>Prediction based anomaly detectionThis is a slightly modified clone of the Kyle Hundman's github repo [telemanom](https://github.com/khundman/telemanom).See also the article [Detecting Spacecraft Anomalies Using LSTMs and Nonparametric Dynamic Thresholding](https://www.groundai.com/project/detecting-spacecraft-anomalies-using-lstms-and-nonparametric-dynamic-thresholding/bib.bib14) in groundai.com.The approach has two components, an LSTM-NN for prediction and a sophisticated thresholding approach to actually detect anomalies. (The authors compared their approach with Gaussian Tail approaches such as z-score in one dimension or elliptic envelopes for multiple dimensions and showed their approach yields better results.)The LSTM part is quite similar to the approach chosen for [LSTM Neural Networks for Anomaly Detection](https://medium.com/datadriveninvestor/lstm-neural-networks-for-anomaly-detection-4328cb9b6e27)<jupyter_code>import os
import logging
import threading
import pandas as pd
from pandas_profiling import ProfileReport

from dqlearn.information.information import data_info
from dqlearn.ui.minimal import print_data_info

import math
import dill
import numpy as np  
import matplotlib.pyplot as plt
import pydot
import seaborn as seabornInstance
import networkx
from sqlalchemy import Column, Integer, String, Float, DateTime, Boolean, func
from iotfunctions import base
from iotfunctions import anomaly
from iotfunctions import bif
from iotfunctions import entity
from iotfunctions import metadata
from iotfunctions.metadata import EntityType
from iotfunctions.db import Database
from iotfunctions.enginelog import EngineLogging
from iotfunctions import estimator, ui, base, bif
from iotfunctions.base import BaseTransformer
from iotfunctions.ui import (UISingle, UIMultiItem, UIFunctionOutSingle,
                 UISingleItem, UIFunctionOutMulti, UIMulti, UIExpression,
                 UIText, UIStatusFlag, UIParameters)
from iotfunctions.enginelog import EngineLogging
from iotfunctions import pipeline as pp
from iotfunctions.stages import DataWriterSqlAlchemy, DataWriterFile
from iotfunctions.pipeline import JobController, DataAggregator

import datetime as dt
from scipy.linalg import norm
import scipy as sp
from scipy import signal
from scipy import linalg
from pyod.models.cblof import CBLOF
from pyod.models.ocsvm import OCSVM
from pyod.models.xgbod import XGBOD

#import pmdarima as pm
import inverse_covariance as icov # from package skggm
from inverse_covariance import QuicGraphLasso

from sklearn.model_selection import train_test_split 
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_squared_error
from sklearn import metrics
from sklearn.decomposition import PCA
from sklearn.decomposition import KernelPCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis as QDA

from skimage import util as skiutil # for nifty windowing

from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error


from keras.models import Sequential
from keras import optimizers, losses, activations, models
from keras.callbacks import ModelCheckpoint, EarlyStopping, LearningRateScheduler, ReduceLROnPlateau
from keras.layers import LSTM, Dense, Flatten, Input, Dropout, Convolution1D, MaxPooling1D, GlobalMaxPooling1D, GlobalAveragePooling1D, concatenate
from keras.utils import plot_model
from sklearn.metrics import accuracy_score, f1_score
from sklearn.model_selection import train_test_split

from IPython.display import SVG
from keras.utils import model_to_dot
import telemanom
from telemanom.helpers import Config
from telemanom.errors import Errors
import telemanom.helpers as helpers
from telemanom.channel import Channel
from telemanom.modeling import Model

%matplotlib inline
# load original data set
#!curl -O https://s3-us-west-2.amazonaws.com/telemanom/data.zip && unzip data.zip -d telemanom && rm data.zip
!mkdir -p ./telemanom/data/2020-01-29_15.00.10/data
# set up the default model parameters
# hidden layers 	2
# units in hidden layers 	80
# sequence length (ls) 	250
# training iterations 	30
# dropout 	0.3
# batch size 	64
# optimizer 	Adam
logger = helpers.setup_logging()
logger.setLevel(logging.DEBUG)
conf = Config("/home/markus/src/mmfunctions/telemanom/config.yaml")
# Load data from 
device="E-3"
chan = Channel(conf, device)
helpers.make_dirs(conf.use_id, conf, "./telemanom")
chan.load_data("./telemanom")
# chan.train
df2 = pd.DataFrame(chan.train)
df2
df2.describe()
# split training data into windows of size 260: 
#   lookback aka l_s (defined in config.yaml as 250) + n_predictions (defined in config.yaml as 10)
#   to predict the last 10 data points from the past 250 ones
#   this is quite similar to our windowsize parameter !
chan.shape_data(chan.train, train=True)
print(chan)
# T-9,MSL,"[[780, 810], [890, 970]]","[point, point]",1096
model = Model(conf, conf.use_id, chan, "./telemanom", False)
model.train_new(chan)
# Stacked LSTM with a Dense Layer of length 10 to predict the next 10 values
print (model)
df = pd.DataFrame(chan.X_train[0])
df.describe()
# P-10,MSL,"[[4590, 4720]]",[point],6100
#df.iloc[2750:2800]
model.save("./telemanom")
# print (chan.X_train)
model.batch_predict(chan, Path="./telemanom")
# predict test data batch-wise and append the results
print (chan.y_hat.shape, chan.y_test.shape)
# smooth the prediction error and apply exponential weights to it
errors = Errors(chan, conf, conf.use_id, "./telemanom")

#  for each overlapping window establish a threshold so that removing error points above it 
# maximizes the reduction of mean and standard deviation. Sort of an adaptive z-score 
errors.process_batches(chan)
print (errors.E_seq, " \n ", errors.anom_scores)<jupyter_output>[(5570, 5709)]  
  [{'start_idx': 5320, 'end_idx': 5389, 'score': 3.504268339986817}, {'start_idx': 5390, 'end_idx': 5459, 'score': 3.7353472227877074}]<jupyter_text>Real life dataNow we try with real life data, just copy the pump data from the internal github into the telemanom subfolder and run the next steps in sequence.We attempt to spot anomalies in `raw_vibrations_x0`, incidentally the first non-array, non-timestamp column of the pump dataset. I'm simply stripping away timestamps, arrays and any vibration related column but the first.<jupyter_code>df_good = pd.read_csv('telemanom/normal_pumps.csv', parse_dates=['RCV_TIMESTAMP_UTC'])
df_bad = pd.read_csv('telemanom/fault_pumps.csv', parse_dates=['RCV_TIMESTAMP_UTC'])
df_ugly = pd.read_csv('telemanom/test_pumps.csv', parse_dates=['RCV_TIMESTAMP_UTC'])

print (df_good.shape, df_bad.shape, df_ugly.shape)

# full list of column names for further reference (aka copy-and-paste)
list_attr = ['raw_vibrations_x_0',
       'raw_vibrations_x_1', 'raw_vibrations_x_2', 'raw_vibrations_x_3',
       'raw_vibrations_x_4', 'raw_vibrations_x_5', 'raw_vibrations_x_6',
       'raw_vibrations_x_7', 'raw_vibrations_x_8', 'raw_vibrations_x_9',
       'raw_vibrations_x_10', 'raw_vibrations_x_11', 'raw_vibrations_x_12',
       'raw_vibrations_x_13', 'raw_vibrations_x_14', 'raw_vibrations_y_0',
       'raw_vibrations_y_1', 'raw_vibrations_y_2', 'raw_vibrations_y_3',
       'raw_vibrations_y_4', 'raw_vibrations_y_5', 'raw_vibrations_y_6',
       'raw_vibrations_y_7', 'raw_vibrations_y_8', 'raw_vibrations_y_9',
       'raw_vibrations_y_10', 'raw_vibrations_y_11', 'raw_vibrations_y_12',
       'raw_vibrations_y_13', 'raw_vibrations_y_14', 'raw_vibrations_z_0',
       'raw_vibrations_z_1', 'raw_vibrations_z_2', 'raw_vibrations_z_3',
       'raw_vibrations_z_4', 'raw_vibrations_z_5', 'raw_vibrations_z_6',
       'raw_vibrations_z_7', 'raw_vibrations_z_8', 'raw_vibrations_z_9',
       'raw_vibrations_z_10', 'raw_vibrations_z_11', 'raw_vibrations_z_12',
       'raw_vibrations_z_13', 'raw_vibrations_z_14', 'rms_x_0', 'rms_x_1',
       'rms_x_2', 'rms_x_3', 'rms_x_4', 'rms_y_0', 'rms_y_1', 'rms_y_2',
       'rms_y_3', 'rms_y_4', 'rms_z_0', 'rms_z_1', 'rms_z_2', 'rms_z_3',
       'rms_z_4', 'accel_speed_0', 'accel_speed_1', 'accel_speed_2',
       'accel_speed_3', 'accel_speed_4', 'accel_power_0', 'accel_power_1',
       'accel_power_2', 'accel_power_3', 'accel_power_4']

# predict column 0, raw vibration x-axis from acceleration speed and power
list_attr = ['raw_vibrations_x_0',
       'accel_speed_0', 'accel_speed_1', 'accel_speed_2',
       'accel_speed_3', 'accel_speed_4', 'accel_power_0', 'accel_power_1',
       'accel_power_2', 'accel_power_3', 'accel_power_4']
df_good.head(2)
infoDF = data_info(df_good)
print_data_info(infoDF)
from dqlearn.ui.minimal import validation_transactional, print_validation
import matplotlib.pyplot as plt

valid_list = validation_transactional(df_good)
print_validation(valid_list, 'Transactional check on sensor data')
plt.show()
dfg = df_good.drop(['Unnamed: 0'],axis=1)
dfg
from dqlearn.ui.minimal import validation_sensor

date_column = 'RCV_TIMESTAMP_UTC'
asset_id_column = 'DEVICE_ID'
sensor_columns = ['raw_vibrations_x_0']
valid_list = validation_sensor(dfg, date_column, asset_id_column, sensor_columns)
print_validation(valid_list, 'FPA Check on Sensor Data')
#profile = ProfileReport(df_good, title='Pandas Profiling Report', html={'style':{'full_width':True}})
#profile.to_notebook_iframe()
np.save("telemanom/data/train/Armstrong.npy", df_good[list_attr].values)
np.save("telemanom/data/Armstrong/bad.npy", df_bad[list_attr].values)
np.save("telemanom/data/test/Armstrong.npy", df_ugly[list_attr].values)
# Load data from 
device="Armstrong"
chan = Channel(conf, device)
helpers.make_dirs(conf.use_id, conf, "./telemanom")
chan.load_data("./telemanom")
# chan.train
df2 = pd.DataFrame(chan.train)
df2
# producing overlapping windows of length 260 for lookback (250) and prediction (10)
chan.shape_data(chan.train, train=True)
# init the Keras double stacked LSTM model
model = Model(conf, conf.use_id, chan, "./telemanom", False)
# drink a coffee - training takes almost 20 minutes
model.train_new(chan)<jupyter_output>Train on 27186 samples, validate on 6797 samples
Epoch 1/35
27186/27186 [==============================] - 106s 4ms/step - loss: 0.0147 - val_loss: 0.0118
Epoch 2/35
27186/27186 [==============================] - 107s 4ms/step - loss: 0.0122 - val_loss: 0.0118
Epoch 3/35
27186/27186 [==============================] - 108s 4ms/step - loss: 0.0121 - val_loss: 0.0118
Epoch 4/35
27186/27186 [==============================] - 109s 4ms/step - loss: 0.0120 - val_loss: 0.0118
Epoch 5/35
27186/27186 [==============================] - 109s 4ms/step - loss: 0.0120 - val_loss: 0.0118
Epoch 6/35
27186/27186 [==============================] - 109s 4ms/step - loss: 0.0120 - val_loss: 0.0118
Epoch 7/35
27186/27186 [==============================] - 111s 4ms/step - loss: 0.0120 - val_loss: 0.0118
Epoch 8/35
27186/27186 [==============================] - 112s 4ms/step - loss: 0.0120 - val_loss: 0.0118
Epoch 9/35
27186/27186 [==============================] - 111s 4ms/step - loss: 0.0120 - val_loss: 0.01[...]<jupyter_text>Training parameters```loss_metric: 'mse'     minimize mean square erroroptimizer: 'adam'      sort of adaptive stochastic gradient descentvalidation_split: 0.2  20% of the data is used for validating (val_loss)dropout: 0.3           ditch 30% of the LSTMs results when minimizing the loss function to avoid overfittinglstm_batch_size: 64    number of training data batches to evaluate per optimizer run to update the model’s parameterspatience: 10           try at least 10 times to decrease val_loss smaller by ...min_delta: 0.0003      ... at least min_delta, else stop, so we get at least 'patience' epochsepochs: 35             no more than 35 passes through the entier training dataset.l_s: 250               lookback: num previous timesteps provided to model to predict future valuesn_predictions: 10      number of steps ahead to predict```This is defined in `telemanom/config.yaml`<jupyter_code>print (model)
# Plot training & validation loss values

# predicting takes more than 40 secs
model.batch_predict(chan, Path="./telemanom")
print (chan.y_hat.shape, chan.y_test.shape)
# smooth the prediction error and apply exponential weights to it
errors = Errors(chan, conf, conf.use_id, "./telemanom")

#  for each overlapping window establish a threshold so that removing error points above it 
# maximizes the reduction of mean and standard deviation. Sort of an adaptive z-score 
errors.process_batches(chan)
print (errors.E_seq, " \n ", errors.anom_scores)
model.save("./telemanom")
# we're doing not too bad it seems

model.model.evaluate(chan.X_test, chan.y_test)
npic = 1
fig, ax = plt.subplots(npic+1, 1, figsize=(20, (npic+1) * 4))
ax[npic-1].set_title('Anomaly detection - forecasting vibration based on acceleration')
ax[npic-1].plot(chan.y_test[:,0], lw=1, color='blue')

#ax.scatter(x_axis, temp_sal_high, lw=8, color='red')
for asc in errors.anom_scores:
    x_axis = np.arange(asc['start_idx'],asc['end_idx'],1)
    y_axis = np.zeros(asc['end_idx'] - asc['start_idx'])
    ax[npic-1].scatter(x_axis,y_axis, lw=5, color='red', zorder=10)
    
ax[npic].set_xlabel('Epoch')
ax[npic].set_ylabel('Loss')
ax[npic].set_title('Model training history')
ax[npic].plot(model.history.history['loss'])
ax[npic].plot(model.history.history['val_loss'],color='green')
ax[npic].legend(['Train', 'Validation'], loc='upper left')

#ax[1].plot(chan.test[:,0], lw=1, color='blue')
#ax[1].set_xlabel('Date')
#ax[1].set_ylabel('Compare with raw training data')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11955805955952294, 'fraction_numerical': 0.04221848247603717, 'mean_word_length': 5.664066309117504, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '324172', 'license': 'MIT', 'n_tokens_mistral': 5057, 'n_tokens_neox': 4615, 'n_words': 1235, 'path': 'Telemanom.ipynb', 'repo_name': 'ankit-jha/mmfunctions', 'size': '13667'}