{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":66052,"databundleVersionId":7284431,"sourceType":"competition"}],"dockerImageVersionId":30615,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## KNN and Isolation Forest models\n\n**Based on the work of:**\n\nhttps://www.kaggle.com/code/victorambonati/unsupervised-anomaly-detection\n\nhttps://www.kaggle.com/code/frankmollard/a-story-about-unsupervised-learning","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport random\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.random.seed(123)\n\n# from sklearn.decomposition import PCA\n# from sklearn.manifold import TSNE\n# import umap\nimport statsmodels.api as sm\n\nfrom sklearn import preprocessing\nfrom sklearn.decomposition import PCA\nfrom sklearn.cluster import KMeans\nfrom sklearn.covariance import EllipticEnvelope\n#from pyemma import msm # not available on Kaggle Kernel\nfrom sklearn.ensemble import IsolationForest\n# from sklearn.svm import OneClassSVM\n\n# import keras as ks\n# import tensorflow as tf\n\nfrom scipy import interpolate\n\n# !pip install trimap==1.0.14\n# import trimap\n\n# !pip install pacmap==0.4\n# import pacmap\n\nimport matplotlib.pyplot as plt \nimport seaborn as sns\nimport plotly.express as px\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"_uuid":"7f97315e-5b74-41b9-bf9b-9d69193973c1","_cell_guid":"75b3086b-fe5e-4008-b104-a0dd893c3e4c","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Read the Data","metadata":{"_uuid":"9b88413d-d07d-4eb4-b8cb-98f3456f8cfc","_cell_guid":"2832914c-442b-406e-9cb7-26702c17c890","trusted":true}},{"cell_type":"code","source":"raw_train_data = pd.read_csv('/kaggle/input/ai-serving-grid-stability/train.csv', index_col='id', parse_dates=['Datum_Uhrzeit_CET'])\nraw_test_data = pd.read_csv('/kaggle/input/ai-serving-grid-stability/test.csv', index_col='id', parse_dates=['Datum_Uhrzeit_CET'])","metadata":{"_uuid":"4857aec2-a60a-430e-89d2-bd6fe0c31e46","_cell_guid":"0ebfde13-e1b5-4714-8337-a75ffad4a2b8","collapsed":false,"execution":{"iopub.status.busy":"2023-12-13T20:22:29.762809Z","iopub.execute_input":"2023-12-13T20:22:29.763404Z","iopub.status.idle":"2023-12-13T20:23:06.121476Z","shell.execute_reply.started":"2023-12-13T20:22:29.763372Z","shell.execute_reply":"2023-12-13T20:23:06.120625Z"},"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"raw_train_data.head(10)","metadata":{"_uuid":"ae9b3c42-295f-4fa9-aaed-5132ac3d919a","_cell_guid":"98fb8541-f1f7-4d73-be04-b23adaf43448","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"raw_test_data.head(10)","metadata":{"_uuid":"412fab6f-6852-4fcb-85c4-c39781e09876","_cell_guid":"f0da6738-7749-49a0-854c-fff95041e2e4","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Explore the Data","metadata":{"_uuid":"b3b45733-69a6-4abe-adca-6202347244ed","_cell_guid":"d90af987-7d7c-432e-ab63-77c3ed4e4fb9","trusted":true}},{"cell_type":"code","source":"# The Dataset contains data from two different control areas. Make sure to separate them when handling the data ;-)\ncontrol_area_1 = raw_train_data[raw_train_data['controlArea']==1]\ncontrol_area_2 = raw_train_data[raw_train_data['controlArea']==2]\n\ncontrol_area_1[['correctedDemand', 'aFRRactivation', 'aFRRrequest']][1000:2000].plot()\ncontrol_area_2[['correctedDemand', 'aFRRactivation', 'aFRRrequest']][1000:2000].plot()","metadata":{"_uuid":"a3770f01-5aed-46aa-8b23-cc8e6aed628e","_cell_guid":"7de9bea4-020a-4743-b720-2232d0ae5201","collapsed":false,"execution":{"iopub.status.busy":"2023-12-13T20:23:06.122689Z","iopub.execute_input":"2023-12-13T20:23:06.123064Z","iopub.status.idle":"2023-12-13T20:23:07.504437Z","shell.execute_reply.started":"2023-12-13T20:23:06.123031Z","shell.execute_reply":"2023-12-13T20:23:07.503296Z"},"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# The test dataset is split into multiple data segments\nraw_test_data['test_data_segment_id'].unique()","metadata":{"_uuid":"61864e8e-35ea-4a8f-8fbd-4715e8f68639","_cell_guid":"399ddd0f-2e4b-4dc6-805a-b7da15fc6e6d","collapsed":false,"execution":{"iopub.status.busy":"2023-12-13T20:23:07.507316Z","iopub.execute_input":"2023-12-13T20:23:07.507664Z","iopub.status.idle":"2023-12-13T20:23:07.517051Z","shell.execute_reply.started":"2023-12-13T20:23:07.507631Z","shell.execute_reply":"2023-12-13T20:23:07.515879Z"},"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# All test data segments have the same length\nraw_test_data.groupby('test_data_segment_id').size()","metadata":{"_uuid":"68b159b9-de13-443f-a3d6-0518009ca053","_cell_guid":"19db0301-008d-414b-862b-38638044a489","collapsed":false,"execution":{"iopub.status.busy":"2023-12-13T20:23:07.518421Z","iopub.execute_input":"2023-12-13T20:23:07.518873Z","iopub.status.idle":"2023-12-13T20:23:07.533792Z","shell.execute_reply.started":"2023-12-13T20:23:07.518832Z","shell.execute_reply":"2023-12-13T20:23:07.532614Z"},"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Each Data Segment belongs to one of the both control Areas\nraw_test_data.groupby('test_data_segment_id').mean()['controlArea']","metadata":{"_uuid":"838ca755-1843-479e-9ca3-ebf5668f22f4","_cell_guid":"bc6789b0-fc76-42ba-b60e-663e093cea9f","collapsed":false,"execution":{"iopub.status.busy":"2023-12-13T20:23:07.535059Z","iopub.execute_input":"2023-12-13T20:23:07.536222Z","iopub.status.idle":"2023-12-13T20:23:07.556977Z","shell.execute_reply.started":"2023-12-13T20:23:07.536156Z","shell.execute_reply":"2023-12-13T20:23:07.555855Z"},"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some features\n\n**Based on the work of:**\n\nhttps://www.kaggle.com/code/victorambonati/unsupervised-anomaly-detection\n\nhttps://www.kaggle.com/code/frankmollard/a-story-about-unsupervised-learning","metadata":{"_uuid":"88104ca0-d79d-4609-a839-18f145924348","_cell_guid":"91ea434d-5af7-4864-8ca2-0998a77b91be","trusted":true}},{"cell_type":"code","source":"# the hours and if it's night or day (6:00-18:00)\ndef feat_eng(df):\n    df['hour'] = df['Datum_Uhrzeit_CET'].dt.hour\n    df['daylight'] = ((df['hour'] >= 6) & (df['hour'] <= 18)).astype(int)\n    # the day of the week (Monday=0, Sunday=6) and if it's a week end day or week day.\n    df['dayofweek'] = df['Datum_Uhrzeit_CET'].dt.dayofweek\n    df['workday'] = (df['dayofweek'] < 5).astype(int)\n    return df","metadata":{"_uuid":"9ececcfd-b9f6-447e-b626-2d8aad6b1751","_cell_guid":"3c6e8a1f-d15b-4e47-8ebb-fbeb21dfe07c","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_1 = feat_eng(control_area_1)\ntrain_2 = feat_eng(control_area_2)\ntest = feat_eng(raw_test_data)\n#\ntest_1 = test[test['controlArea']==1]\ntest_2 = test[test['controlArea']==2]","metadata":{"_uuid":"faa818d3-4c9f-4d34-9053-0fde3561c94a","_cell_guid":"0da5736d-8531-4ef8-9c32-5a6d6b615f69","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_1.describe()","metadata":{"_uuid":"2a1d0d0a-2516-46fc-947d-4822f334b445","_cell_guid":"a77734d1-b729-4731-ad99-a18c0de1d353","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_2.describe()","metadata":{"_uuid":"34705249-a5e7-43d8-b874-0ae317ed9e6f","_cell_guid":"2ddf18ae-a1e4-44c5-ba86-903dcfda46c3","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_1.describe()","metadata":{"_uuid":"6c5e12a2-173c-4804-83ab-d0fdbda716f3","_cell_guid":"236f3eff-f086-4152-b557-c70cf96de109","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# KNN\n\nIf the point is to far from the centroid of its group, the point is an outlier.","metadata":{"_uuid":"42d9b217-3bb4-4628-8968-4965098a5351","_cell_guid":"84391e03-4dea-4984-a13a-381d73afb033","trusted":true}},{"cell_type":"code","source":"# return Series of distance between each point and his distance with the closest centroid\ndef getDistanceByPoint(data, model):\n    distance = pd.Series()\n    for i in range(0,len(data)):\n        Xa = np.array(data.loc[i])\n        Xb = model.cluster_centers_[model.labels_[i]-1]\n        distance.loc[i] = np.linalg.norm(Xa-Xb)\n    return distance\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Control Area 1\nX = train_1.drop(columns=[\"Datum_Uhrzeit_CET\", \"controlArea\"])\nss_scaler_1 = preprocessing.StandardScaler()\nnp_scaled = ss_scaler_1.fit_transform(X)\ndata_1 = pd.DataFrame(np_scaled)\n\n# Control Area 2\nX = train_2.drop(columns=[\"Datum_Uhrzeit_CET\", \"controlArea\"])\nss_scaler_2 = preprocessing.StandardScaler()\nnp_scaled = ss_scaler_2.fit_transform(X)\ndata_2 = pd.DataFrame(np_scaled)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# It takes hours...\n# calculate with different number of centroids to see the loss plot (elbow method)\n# n_cluster = range(5, 15)\n# kmeans = [KMeans(n_clusters=i).fit(data) for i in n_cluster]\n# scores = [kmeans[i].score(data) for i in range(len(kmeans))]\n# fig, ax = plt.subplots()\n# ax.plot(n_cluster, scores)\n# plt.show()\n\n# Set number of groups equal to 10. \nnn = 10\n\nkmeans_1 = KMeans(n_clusters=nn).fit(data_1)\nkmeans_2 = KMeans(n_clusters=nn).fit(data_2)\n\ncluster = pd.DataFrame()\ncluster['cluster'] = kmeans_1.predict(data_1)\ncluster['cluster'].value_counts()\n\n# TEST 1\nX = test_1.drop(columns=[\"Datum_Uhrzeit_CET\", \"controlArea\", \"test_data_segment_id\"])\nnp_scaled = ss_scaler_1.transform(X)\ndata = pd.DataFrame(np_scaled)\ncluster = pd.DataFrame()\ncluster['cluster'] = kmeans_1.predict(data)\ncluster['cluster'].value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get the distance between each point and its nearest centroid. The biggest distances are considered as anomaly\ndistance = getDistanceByPoint(data, kmeans_1)\ndistance.describe()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define the proportion of outliers\n# nn = 10\n# outliers_fraction = 0.1 # 0.2538\n# outliers_fraction = 0.2 # 0.3917\n# outliers_fraction = 0.3 # 0.4442\noutliers_fraction = 0.5 # 0.5157\n# outliers_fraction = 0.6 # 0.5117\n\nnumber_of_outliers = int(outliers_fraction*len(distance))\nthreshold = distance.nlargest(number_of_outliers).min()\n# contain the anomaly result of method 2.1 Cluster (0:normal, 1:anomaly) \ndata['anomalyKNN'] = (distance >= threshold).astype(int)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def knn_outlier(test, ss_scaler, kmeans, outliers_fraction):\n    X = test.drop(columns=[\"Datum_Uhrzeit_CET\", \"controlArea\", \"test_data_segment_id\"])\n    np_scaled = ss_scaler.transform(X)\n    data = pd.DataFrame(np_scaled)\n    # get the distance between each point and its nearest centroid. The biggest distances are considered as anomaly\n    distance = getDistanceByPoint(data, kmeans)\n    number_of_outliers = int(outliers_fraction*len(distance))\n    threshold = distance.nlargest(number_of_outliers).min()\n    data['outlier'] = (distance >= threshold).astype(int)\n    return data\n\n","metadata":{"_uuid":"155be6a1-424b-49b9-8fc0-fc98a2e4c9ef","_cell_guid":"58fdbdf5-28e5-4134-8cc7-fd21d8989ad5","collapsed":false,"jupyter":{"outputs_hidden":false},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Isolation Forest\n\nIsolation Forest is an algorithm for data anomaly detection initially developed by Fei Tony Liu in 2008. Isolation Forest detects anomalies using binary trees. The algorithm has a linear time complexity and a low memory requirement, which works well with high-volume data. (https://en.wikipedia.org/wiki/Isolation_forest)","metadata":{}},{"cell_type":"code","source":"# Control Area 1\nX = train_1.drop(columns=[\"Datum_Uhrzeit_CET\", \"controlArea\"])\nss_scaler_1 = preprocessing.StandardScaler()\nnp_scaled = ss_scaler_1.fit_transform(X)\ndata_1 = pd.DataFrame(np_scaled)\n\n# Control Area 2\nX = train_2.drop(columns=[\"Datum_Uhrzeit_CET\", \"controlArea\"])\nss_scaler_2 = preprocessing.StandardScaler()\nnp_scaled = ss_scaler_2.fit_transform(X)\ndata_2 = pd.DataFrame(np_scaled)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"outliers_fraction = .5\nmodel_1 =  IsolationForest(contamination = outliers_fraction)\nmodel_1.fit(data_1)\n# add the data to the main  \npred = pd.DataFrame()\npred['anomaly'] = pd.Series(model_1.predict(data_1))\npred['anomaly'] = pred['anomaly'].map( {1: 0, -1: 1} )\nprint(pred['anomaly'].value_counts())\n\nmodel_2 =  IsolationForest(contamination = outliers_fraction)\nmodel_2.fit(data_2)\n# add the data to the main  \npred = pd.DataFrame()\npred['anomaly'] = pd.Series(model_2.predict(data_2))\npred['anomaly'] = pred['anomaly'].map( {1: 0, -1: 1} )\nprint(pred['anomaly'].value_counts())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TEST 1\nX = test_1.drop(columns=[\"Datum_Uhrzeit_CET\", \"controlArea\", \"test_data_segment_id\"])\nnp_scaled = ss_scaler_1.transform(X)\ndata = pd.DataFrame(np_scaled)\npred = pd.DataFrame()\npred['anomaly'] = pd.Series(model_1.predict(data))\npred['anomaly'] = pred['anomaly'].map( {1: 0, -1: 1} )\nprint(pred['anomaly'].value_counts())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef isof(test, ss_scaler, model):\n    X = test.drop(columns=[\"Datum_Uhrzeit_CET\", \"controlArea\", \"test_data_segment_id\"])\n    np_scaled = ss_scaler.transform(X)\n    data = pd.DataFrame(np_scaled)\n    data['anomaly'] = pd.Series(model.predict(data))\n    data['anomaly'] = data['anomaly'].map( {1: 0, -1: 1} )\n    return data","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission KNN","metadata":{}},{"cell_type":"code","source":"# Submission?\nSUB = True\n\nif SUB:\n    dataset_ids = raw_test_data['test_data_segment_id'].unique()\n\n    submission = np.array([])\n\n    for dataset_id in dataset_ids:\n        sample = raw_test_data[raw_test_data['test_data_segment_id'] == dataset_id]\n        if np.any(sample.controlArea == 1):  \n            data_1 = knn_outlier(sample, ss_scaler_1, kmeans_1, outliers_fraction)\n            submission = np.append(submission, data_1['outlier'])\n        else: \n            data_2 = knn_outlier(sample, ss_scaler_2, kmeans_2, outliers_fraction)\n            submission = np.append(submission, data_2['outlier'])\n            \n    submission = pd.DataFrame({'anomaly': submission})\n    submission.index.name = 'id'\n    # save as submission.csv in order to submit via this kaggle notebook\n    print(submission['anomaly'].sum()/submission.shape[0])\n    submission.to_csv('submission_knn.csv')\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission Isol-Forest","metadata":{}},{"cell_type":"code","source":"SUB = True\n\nif SUB:\n    dataset_ids = raw_test_data['test_data_segment_id'].unique()\n\n    submission = np.array([])\n\n    for dataset_id in dataset_ids:\n        sample = raw_test_data[raw_test_data['test_data_segment_id'] == dataset_id]\n        if np.any(sample.controlArea == 1):  \n            data_1 = isof(sample, ss_scaler_1, model_1)\n            submission = np.append(submission, data_1['anomaly'])\n        else: \n            data_2 = isof(sample, ss_scaler_2, model_2)\n            submission = np.append(submission, data_2['anomaly'])\n            \n    submission = pd.DataFrame({'anomaly': submission})\n    submission.index.name = 'id'\n    # save as submission.csv in order to submit via this kaggle notebook\n    print(submission['anomaly'].sum()/submission.shape[0])\n    submission.to_csv('submission_isol.csv')\n\n","metadata":{},"execution_count":null,"outputs":[]}]}