Isolation forest + lstm autoencoder anomaly detection

1. isolation forest 로 정상에 가까운 클러스터, 비정상에 가까운 클러스터 이렇게 두개 묶음 만들기

2.  데이터 train_정상/train_비정상/test_정상/test_비정상 이케 데이터 나눠준다

3. train_정상 만 가지고 lstm autoencoder 학습시키고 테스트 데이터 넣어서 threshold 를 구해준다 

=> 이 스레숄드가 뭐나면? 재구성손실값임. 한마디로 정상치만 가지고 학습시켰을때 똑같이 정상치라면 재구성할떄 오류가 많이 나지않음.스레숄드가 낮게나온다는 뜻임. 그렇다면 정상, 그런데 완전 다른 그림을 학습시킨다면 재구성할때 오차가 엄청날것임. 스레숄드가 엄청 크게 나온다는뜻. 스레숄드보다 높은값이 나온것은 정상이랑 많이 다른양상이라는 뜻이므로 anomaly 로 구분한다.

4. 이렇게 전체 데이터를 집어넣어봐서 스레숄드보다 높은값을 보이면 (재구성손실이 많다면) anomaly 로 라벨링하고, 스레숄드보다 낮은값을 보이면(재구성손실이 적다면, 정상치 트레인데이터와 비슷한양상을 보인다면) 정상값으로 라벨링한다.

스레숄드 구하는법은 아래와같음

** 정밀도 ( precision) : 참이라고 예측한 데이터중 실제로 참인 데이터 

**제현율 ( RECALL ) : 실제로 참인 데이터중에서 참이라고 예측한 데이터

긁어옴)) : 재현율이 높지만 정밀도가 낮다는 것은 찾은 데이터의 수가 많지만 그 중 실제 찾으려는 대상의 비율은 낮았다는 것을 의미합니다. 이와 반대로 정밀도가 높지만 재현율이 낮다는 것은 찾은 데이터의 수는 작지만 그 중 실제 찾으려는 대상의 비율은 높다는 것을 의미합니다. 이상적으로는, 정밀도와 재현율이 모두 높은 것이 좋습니다. 이는 찾은 데이터 수가 많으며 그 중 실제 찾으려는 대상의 비율이 높다는 것을 의미합니다. 그렇기 때문에 threshold 값은 이 정밀도와 재현률이 둘다 좋은지점을 찾아 threshold 라고 정의해줍니다.

 

 

참조 : 

https://www.koreascience.or.kr/article/CFKO202125036427364.pdf (국내학회)

https://towardsdatascience.com/lstm-autoencoder-for-extreme-rare-event-classification-in-keras-ce209a224cfb

https://velog.io/@jaehyeong/LSTM-Autoencoder-for-Anomaly-Detection 

 

 

 

 

 

 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import math
from pylab import rcParams
from collections import Counter
import tensorflow as tf
from tensorflow.keras import Model ,models, layers, optimizers, regularizers
from tensorflow.keras.callbacks import ModelCheckpoint
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn import metrics
%matplotlib inline
from sklearn.ensemble import IsolationForest
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from mpl_toolkits.mplot3d import Axes3D
import pickle
import glob
from tensorflow.python.keras.metrics import Metric
import sklearn

# lsmt ae 에서 사용할 함수 
# LSTM 모델은 (samples, timesteps, feature)에 해당하는 3d 차원의 shape을 가지므로, 데이터를 시퀀스 형태로 변환한다.
def temporalize(X, y, timesteps): 
	output_X = []  
	output_y = []
	for i in range(len(X) - timesteps - 1):  
		t = []
		for j in range(1, timesteps + 1):
			# Gather the past records upto the lookback period
			t.append(X[[(i + j + 1)], :])
		output_X.append(t)
		output_y.append(y[i + timesteps + 1])
	return np.squeeze(np.array(output_X)), np.array(output_y)

content={}
#데이터별로 이상치 index를 담아줄 딕셔너리. 마지막에 anomaly라고 라벨링된 index만 뽑아서 "데이터이름" : [1,2,4...] 이런식으로 key,value 값을 저장해줄것.
parameter = [0.01 , 0.015, 0.02, 0.025, 0.03]

for k in filenames:
    for j in parameter:
        print(j)
        df = pd.read_csv(str(k))
        
        
        # 1차 라벨링 부분
        #isolationForest 적용 + 그림으로 보이기 ---------------------------------------------------------
        
        from sklearn.ensemble import IsolationForest
        clf=IsolationForest(n_estimators=50, max_samples=50, contamination = float(j), 
                                max_features=1.0, bootstrap=False, n_jobs=-1, random_state=None, verbose=0)
        clf.fit(df)
        pred = clf.predict(df)
        df['anomaly'] = pred  # isolation forest로 1차 라벨링된것을 새피처로써 추가함.
        outliers = df.loc[df['anomaly'] == -1]
        outlier_index = list(outliers.index) # 1차 라벨링해서 이상치값의 index 만 outlier_index 리스트로 만듬. 
        pca = PCA(n_components = 5) # isolation 적용한후, 피처개수가 너무많아서 차원축소해준다. 

        scaler = StandardScaler()
        X = scaler.fit_transform(df)
        X_reduce = pca.fit_transform(X)
        fig = plt.figure()
        ax = fig.add_subplot(111)
        ax.scatter(X_reduce[:, 0], X_reduce[:, 1], s=4, lw=1, label="inliers",c="green")
        ax.scatter(X_reduce[outlier_index,0],X_reduce[outlier_index,1],
                   lw=1, s=4, c="red", label="outliers")
        ax.legend()
        plt.show()
        print("isolation forest 적용후값  : " + str(Counter(pred)))
        
         
        # 2차 라벨링 부분 
        # Train/valid data split -------------------------------------------------------------------------------
        
        input_x = X_reduce   # isolation forest 로 1차 라벨링한값을 차원축소한후 lstm ae 의 인풋데이터로 사용할것
        input_y = []
        # 1차 라벨링에서 1 정상 -1 비정상으로 1차 라벨링된값을 => 1 비정상 0 정상으로 숫자만 바꿔줌
        for i in range(len(input_x)):   # 1 , -1
            if i in outlier_index:
                input_y.append(1) # 비정상
            else :
                input_y.append(0) # 정상
        input_y = np.array(input_y)  
        
        timesteps = 3
        x, y = temporalize(input_x, input_y, timesteps)   # LSTM 모델은 (samples, timesteps, feature)에 해당하는 3d 차원의 shape을 가지므로, 데이터를 시퀀스 형태로 변환한다.


        # Split into train and valid
        x_train, x_valid, y_train, y_valid = train_test_split(x, y, test_size=0.2)
        n_features = input_x.shape[1]
        
        # For training the autoencoder, split 0 / 1
        x_train_y0 = x_train[y_train == 0] # 정상데이터 y0
        x_train_y1 = x_train[y_train == 1] # 비정상데이터 y1
        x_valid_y0 = x_valid[y_valid == 0]
        x_valid_y1 = x_valid[y_valid == 1]

        x_y0 = x[y==0] # 정상데이터 y0
        x_y1 = x[y==1] # 비정상데이터 y1
        
        def flatten(X):
            flattened_X = np.empty((X.shape[0], X.shape[2]))  # sample x features array.
            for i in range(X.shape[0]):
                flattened_X[i] = X[i, (X.shape[1]-1), :]
            return(flattened_X)
        def scale(X, scaler):
            for i in range(X.shape[0]):
                X[i, :, :] = scaler.transform(X[i, :, :])  
            return X

        scaler = StandardScaler().fit(flatten(x_train_y0))
        x_train_y0_scaled = scale(x_train_y0, scaler)
        x_valid_scaled = scale(x_valid, scaler)
        x_valid_y0_scaled = scale(x_valid_y0, scaler)
        x_total_scaled = scale(x, scaler)   
        x_y0_sacled = scale(x_y0, scaler)

        # LSTM Autoencoder 적용 -------------------------------------------------------------------------------
        epochs = 5
        batch = 128
        lr = 0.001
        lstm_ae = models.Sequential()
        # Encoder
        lstm_ae.add(layers.LSTM(32, activation='relu', input_shape=(timesteps, n_features), return_sequences=True))
        lstm_ae.add(layers.LSTM(16, activation='relu', return_sequences=False))
        lstm_ae.add(layers.RepeatVector(timesteps))
        # Decoder
        lstm_ae.add(layers.LSTM(16, activation='relu', return_sequences=True))
        lstm_ae.add(layers.LSTM(32, activation='relu', return_sequences=True))
        lstm_ae.add(layers.TimeDistributed(layers.Dense(n_features)))
        # compile
        lstm_ae.compile(loss='mse', optimizer=optimizers.Adam(lr))
        # fit
        history = lstm_ae.fit(x_y0_sacled, x_y0_sacled,
                             epochs=epochs, batch_size=batch,
                             validation_data=(x_valid_y0_scaled, x_valid_y0_scaled))
        
        
        
        # normal/anomaly 를 구분해주는 기준 threshold 를 구함  -------------------------------------------------------

        valid_x_predictions = lstm_ae.predict(x_valid_scaled)
        mse = np.mean(np.power(flatten(x_valid_scaled) - flatten(valid_x_predictions), 2), axis=1)

        error_df = pd.DataFrame({'Reconstruction_error':mse, 
                                 'True_class':list(y_valid)})
        precision_rt, recall_rt, threshold_rt = sklearn.metrics.precision_recall_curve(error_df['True_class'], error_df['Reconstruction_error'])
        
        
        # precision/recall 값
        index_cnt = [cnt for cnt, (p, r) in enumerate(zip(precision_rt, recall_rt)) if p==r][0]
        #print('precision: ',precision_rt[index_cnt],', recall: ',recall_rt[index_cnt])

        # threshold 값
        threshold_fixed = threshold_rt[index_cnt]
        #print('threshold: ',threshold_fixed)


        # threshold 보다 낮으면 정상(normal), 높으면 이상(anomaly)으로 판단한다.
        # 정밀도 ( precision) : 참이라고 예측한 데이터중 실제로 참인 데이터 (참 / 참인데참,거짓인데참)
        # 제햔율 ( RECALL ) : 실제로 참인 데이터중에서 참이라고 예측한 데이터? (참/ 참인데참, 참인데 거짓))

        
        
        # threshold 값을 기준으로 normal/anomlay 인지 라벨링해줌 ----------------------------------------------------------------
        test_x_predictions = lstm_ae.predict(x_total_scaled) 
        mse = np.mean(np.power(flatten(x_total_scaled) - flatten(test_x_predictions), 2), axis=1)
        error_df = pd.DataFrame({'Reconstruction_error': mse,
                                 'True_class': y.tolist()})
        groups = error_df.groupby('True_class')
        fig, ax = plt.subplots()

        for name, group in groups:
            ax.plot(group.index, group.Reconstruction_error, marker='o', ms=3.5, linestyle='',
                    label= "Break" if name == 1 else "Normal")
        ax.hlines(threshold_fixed, ax.get_xlim()[0], ax.get_xlim()[1], colors="r", zorder=100, label='Threshold')
        ax.legend()
        plt.title("Reconstruction error for different classes")
        plt.ylabel("Reconstruction error")
        plt.xlabel("Data point index")
        plt.show();
        
        
        # classification by threshold
        pred_y = [1 if e > threshold_fixed else 0 for e in error_df['Reconstruction_error'].values]
        LABELS = ['Normal', 'Break']
        conf_matrix = sklearn.metrics.confusion_matrix(error_df['True_class'], pred_y)
        plt.figure(figsize=(7, 7))
        sns.heatmap(conf_matrix, xticklabels=LABELS, yticklabels=LABELS, annot=True, fmt='d')
        plt.title('Confusion Matrix')
        plt.xlabel('Predicted Class'); plt.ylabel('True Class')
        plt.show()


        # TN/FP/FN 구하기 ------------------------------------------------------------------------------------
        # 여기서 true 값은 1차 라벨링값 predict 는 2차 라벨링값
        
        df_p = pd.DataFrame(pred_y, columns=['anomaly'])
        df_t = pd.DataFrame(error_df['True_class'].tolist(), columns=['anomaly'])
        
        anomaly_list_p = df_p[df_p['anomaly']==1].index
        anomaly_list_T = df_t[df_t['anomaly']==1].index
        anomaly_list_p = (anomaly_list_p +2).to_list()
        anomaly_list_T = (anomaly_list_T +2).to_list()

        TN = list(set(anomaly_list_T) & set(anomaly_list_p)) 
        FP = list(set(anomaly_list_p) - set(anomaly_list_T))
        FN = list(set(anomaly_list_T) - set(anomaly_list_p))
        
                          
        # anomaly index 만 리스트로 만들어서 dictionary 에 value 로 저장. ------------------------------------------------------------------------------------------
        isolationforest_parameter = j
        str_isolationforest_parameter = str(isolationforest_parameter)[2:]
        if(len(str_isolationforest_parameter) != 5  ):
            for i in range(5 - len(str_isolationforest_parameter)):
                str_isolationforest_parameter=str_isolationforest_parameter + '0'
            
        
        
        Save_dict_name = str(k)[7:20]
        content[Save_dict_name +'_p_'+ str(str_isolationforest_parameter)] = TN+FP+FN
        
        
        with open('anomaly_index.pickle','wb') as fw:
            pickle.dump(content, fw)
            
        file = open("anomaly_index.pickle",'rb')
        content = pickle.load(file)
        content.keys()