{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":66052,"databundleVersionId":7284431,"sourceType":"competition"},{"sourceId":7613523,"sourceType":"datasetVersion","datasetId":4433676}],"dockerImageVersionId":30646,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Introduction\n\nThe objective of this project is to identify anomalies within the data from the PICASSO energy grid balancing system.\nThe problem is an unsupervised, multi-variate timeseries anomaly detection problem. We have elected to employ a hybrid strategy that combines technically motivated indicators for anomalies (\"rule-based\"-detecton of abnormal system states), along with a purely data-driven machine learning system for identification of statistical anomalies in the mutivariate feature space.\n\nThe reasoning for this strategy is as follows: A machine learning approach will efficiently detect statisticaly anomalies (not only point-like anomalies, but also complex contextual multivariate anomalies), but will likely not account for all abnormal system states, that might be categorized as \"anomalies\" in the test set. For example, sensor transmission errors might be regraded as abnormal system states, but occur frequently, and thus, won't be detected as statistical anomalies by an machine learning approach. Employing a rule-based method, inspired by both the functional relationship between features  as well as \"known\" abnormal states of the system, provides a comprehensive capture of abnormal system states results. In combination both approaches might result in a very good caputure of the defintion of \"anomalies\" in this challenge.\n\n\nBelow, the images illustrate our modeling setup, encompassing several pre- and post-processing routines. These routines aim to both amplify our model's performance and rectify the deficiencies inherent with a solely statistical approach.\n\n![Image](https://i.imgur.com/MdJbKm6.png)\n\nIn the following the individual components of the approach are described in detail.","metadata":{}},{"cell_type":"markdown","source":"# Setup\n## Required packages","metadata":{}},{"cell_type":"code","source":"# System libaries\nimport os\nimport sys\nimport gc  # garbage collector interface\nimport random\nimport time\nimport warnings\n\n# Other libraries\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom IPython.display import display_html\n\nfrom scipy import signal\nfrom sklearn.preprocessing import StandardScaler\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras.layers import Input, LSTM, RepeatVector, TimeDistributed, Dense, Lambda, Dropout\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras import backend as K\nfrom tensorflow.keras import regularizers\nfrom tensorflow.keras.callbacks import Callback\n\n# Pandas configurations\npd.options.mode.chained_assignment = None\npd.set_option('display.max_rows', 500)  # Set maximum number of rows displayed\nwarnings.simplefilter(action='ignore', category=FutureWarning)","metadata":{"execution":{"iopub.status.busy":"2024-02-13T17:34:17.004546Z","iopub.execute_input":"2024-02-13T17:34:17.004954Z","iopub.status.idle":"2024-02-13T17:34:30.756362Z","shell.execute_reply.started":"2024-02-13T17:34:17.004921Z","shell.execute_reply":"2024-02-13T17:34:30.755535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load datasets","metadata":{}},{"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'])\n\nX_train = raw_train_data.copy()\nX_test = raw_test_data.copy()","metadata":{"execution":{"iopub.status.busy":"2024-02-13T17:34:31.070251Z","iopub.execute_input":"2024-02-13T17:34:31.071634Z","iopub.status.idle":"2024-02-13T17:35:06.693904Z","shell.execute_reply.started":"2024-02-13T17:34:31.071589Z","shell.execute_reply":"2024-02-13T17:35:06.692839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Rule-based approach\n\nThrough extensive exploration of the PICASSO working mechanims and explorative analysis of the data provided, we came up with the follwing set of potentially abnormal system stages:\n\n1. Band limit violations for band limited corrected demand\n1. Inconsistencies of demand, correction and corrected demand\n1. Inconsistencies of LFC input + aFRR activation and corrected demand\n1. Extreme values\n1. Extreme movements\n1. (Near) constant values\n1. Inconsistent delays of correction echo\n1. Asynchronous activation of aFRR\n1. Shifted correlation of aFRRactivation and aFRRrequest\n","metadata":{}},{"cell_type":"markdown","source":"## Set parameters for rule-based abnormal state detection","metadata":{}},{"cell_type":"code","source":"# Anomaly type I: Band limit violations for band limited corrected demand\n# Define threshold for consistency of band limits and band limited corrected demand\nband_limit_violation_thr = 0.1\n\n\n# Anomaly type II: Inconsistencies of demand, correction and corrected demand\n# Define threshold for consistency of demand, demand correction, and corrected demand\ndemand_correction_incons_abs_thr = 0.1\n\n\n# Anomaly type III: Inconsistencies of LFC input + aFRR activation and corrected demand\n# Define threshold for consistencies of LFC input + aFRR activation and corrected demand\nLFCInput_incons_thr_quantile_p = 1.0 - 0.001\n\n\n# Anomaly type IV: Extreme values\n# Define thresholds for extremes of general feature values\n# Set p value for quantile of lower / upper value threshold\nextreme_values_lower_thr_quantile_p = 0.001\nextreme_values_upper_thr_quantile_p = 1.0 - extreme_values_lower_thr_quantile_p\n# Directly define thresholds for special cases\nextreme_values_lower_special_thr = pd.Series(data=[np.nan, 1.0, 2.0], index=['controlBandNeg', 'controlBandPos', 'controlBandWidth'])\nextreme_values_upper_special_thr = pd.Series(data=[-1.0, np.nan, np.nan], index=['controlBandNeg', 'controlBandPos', 'controlBandWidth'])\n# Define list of columns to be checked for extreme values\nselect_cols_extreme_values = [ \\\n    'Demand', 'correction', 'correctedDemand', 'FRCE', \\\n    'controlBandPos', 'controlBandNeg', 'controlBandWidth', 'LFCInput', 'aFRRactivation', \\\n    'aFRRrequest', 'correctionEcho', \\\n    'BandLimitedCorrectedDemand']\n\n\n# Anomaly type V: Extreme movements\n# Define threshold for extremes of general feature movements\n# Set p value for quantile of absolute large movement (\"jump\") threshold\nextreme_jumps_abs_thr_quantile_p = 1.0 - 0.002\n# Define list of columns to be checked for extremes\nselect_cols_extreme_jumps = [ \\\n    'Demand', 'correction', 'correctedDemand', 'FRCE', \\\n    'LFCInput', 'aFRRactivation', \\\n    'aFRRrequest', 'correctionEcho']\n\n\n# Anomaly type VI: Near constant values\n# Define threshold value for absolute small movements (i.e. (near) constant values)\nnear_const_abs_thr = 1e-10\n# Define threshold number of adjacent small movements\nnear_const_abs_number_thr = 15\n# Define list of columns to be checked for small movements\nselect_cols_near_const_values = [ \\\n    'Demand', 'correction', 'correctionEcho', 'FRCE', \\\n    'LFCInput', 'aFRRactivation', \\\n    'aFRRrequest']\n# Define list of columns for which ranges of zeros shall not be flagged as anomaly\nselect_cols_near_const_values_elim_0 = [ \\\n    'correction', 'correctionEcho', 'FRCE']\n\n\n# Anomaly type VII: Inconsistent delays of correction echo\n# Define start and end of test shifts for correction time series\n# (anomaly detection for correction and correction echo)\n# to try and match the correction echo\nstart_shift_corr = 0\nend_shift_corr = 5\n# Define threshold for deviation of shifted correction and echo\ncorr_echo_deviation_thr = 1e-3\n# Define threshold for deviation from the most frequent shift with sub threshold deviation\ncorr_echo_deviation_shift_most_frequent_thr = 1\n# Define most frequent shifts as determined from training data\ncorr_echo_deviation_shift_most_frequent = [2, 1]\n# Use training data to calculate most frequent shifts\ncorr_echo_incons_delays_train = False\n\n\n# Anomaly type VIII: Asynchronous activation of aFRR\n# Define thresholds for detection of asynchronous activation \n# Set p value for quantile of absolute activation error threshold (i.e. all aes above the threshold to be marked as anomaly)\nae_thr_quantile_p = 1.0 - 0.005\n# Set distance threshold, i.e. absolute activation errors below this value are considered \"near\"\nae_near_abs_thr = 15.0\n# Set p value for quantile of movement (\"jump\") threshold as basis for the qualification of near values \nae_jump_abs_thr_quantile_p = 1.0 - 0.005\n# Define window size for moving average approach to simulate aFRRactivation\nwindow_size_aFRRactivationSim = [81, 16]\n# Decide if run with training data is to be performed as well\n# (\"True\" leads to potentially huge storage space requirements and long runtime)\nae_include_train = False\n\n\n# Anomaly type IX: Shifted correlation of aFRRactivation and aFRRrequest\nwindow_size_corr = 400\nmin_shift_corr = 0\nmax_shift_corr = 100\ncorr_thr = 0.65\ncorr_noise_thr = 20","metadata":{"execution":{"iopub.status.busy":"2024-02-13T09:53:12.796007Z","iopub.execute_input":"2024-02-13T09:53:12.796307Z","iopub.status.idle":"2024-02-13T09:53:12.808921Z","shell.execute_reply.started":"2024-02-13T09:53:12.796282Z","shell.execute_reply":"2024-02-13T09:53:12.807847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define utility functions","metadata":{}},{"cell_type":"code","source":"def detect_constant_value_ranges(df, near_const_abs_thr, near_const_abs_number_thr, \\\n                                select_cols_near_const_values, select_cols_near_const_values_elim_0=[], \\\n                                fill_first_entry=(near_const_abs_thr + 1.0), flag_first_const_value=False, diff_first_const_value=1):\n    # Assume that the first column differences with the record immediately before the beginning of the time series is above threshold\n    # (Differences are taken with respect to the preceding record, i.e.: a.diff[1] = a[1] - a[0])\n    aux = df.loc[:, select_cols_near_const_values].diff().fillna(fill_first_entry)\n\n    if (flag_first_const_value):\n       aux[aux[select_cols_near_const_values] == 1] = 0\n\n    aux = aux.abs()\n\n    # For each column, mark all first differences below threshold with 1, mark all other differences with 0\n    aux = (aux < near_const_abs_thr).astype(int)\n\n    aux1 = aux.copy()\n    for i in range(1, near_const_abs_number_thr):\n        # For each record, add the marks of the (near_const_abs_number_thr - 1) following records to\n        # obtain a maximum value of near_const_abs_number_thr\n        # The records marked with that maximum value have a sub-threshold difference with their preceding record and\n        # at least the (near_const_abs_number_thr - 1) following records\n        aux1 = aux1.add(aux.shift(-i, fill_value=0))\n\n    # Set all records with the maximum value of near_const_abs_number_thr to 1,\n    # the other records to zero\n    aux1 = (aux1 == near_const_abs_number_thr).astype(int)\n\n    aux = aux1.copy()\n\n    for i in range(1, near_const_abs_number_thr):\n        # Add marks from the (near_const_abs_number_thr - 1) previous records\n        # As a result, all records with a value > 0 are part of a range of nearly constant values\n        # with at least near_const_abs_number_thr differences sub threshold\n        aux1 = aux1.add(aux.shift(i, fill_value=0))\n\n    # For selected columns, prevent ranges of zeros from being flagged as anomalies \n    if (len(select_cols_near_const_values_elim_0) > 0):\n        for column in select_cols_near_const_values_elim_0:\n            aux1.loc[df[column] == 0, column] = 0\n\n    # Detect and mark near constant values\n    # Use OR combination of columns, i.e. a True flag means that\n    # at least one of the columns has a record that is part of near_const_abs_number_thr nearly constant values\n    return((aux1 > 0).astype(int).sum(axis=1).astype(bool))\n\ndef get_window(timeseries1, timeseries2, idx, window_size):\n    window_size = int(window_size/2)\n    window1 = timeseries1[idx-window_size:idx+window_size]\n    window2 = timeseries2[idx-window_size:idx+window_size]\n    return window1, window2\n\ndef get_correlations_for_shifts(\n        window1, \n        window2, \n        min_shift = 0, \n        max_shift = 100\n):\n\n    size = int(len(window1)/2)\n\n    # compute correlation for all shifts\n    c1 = signal.correlate(window2, window1, mode='same')\n    c2 = signal.correlate(window1, window1, mode='same')[size]\n    c3 = signal.correlate(window2, window2, mode='same')[size]\n    corrs = c1 / np.sqrt(c2 * c3)\n\n    # reduce to allowed shifts\n    corrs = corrs[size + min_shift : size + max_shift]\n\n    return corrs\n\ndef get_optimal_shift_and_according_corr(corrs):\n    corr = np.max(corrs)\n    shift = min_shift_corr + np.argmax(corrs)\n    return corr, shift\n\ndef interpolate_nan(data):\n    df = pd.DataFrame({\"data\" : data})\n    df = df.interpolate(\"linear\")\n    df.loc[df.data.isna(), \"data\"] = df.loc[df.data.notna(), \"data\"].to_list()[0]\n    return df.data.to_list()\n\ndef get_shifed_corr_series(\n    timeseries1, \n    timeseries2, \n    window_size = 400,\n    min_shift = 0, \n    max_shift = 100,\n    step = 20,\n    ):\n\n    ts_corr = [None] * len(timeseries1)\n    ts_shift = [None] * len(timeseries1)\n    ts_mean = [None] * len(timeseries1)\n    for idx in range(int(window_size/2), len(timeseries1)-int(window_size/2), step):\n        window1, window2 = get_window(timeseries1, timeseries2, idx, window_size)\n        corrs = get_correlations_for_shifts(window1, window2)\n        corr, shift = get_optimal_shift_and_according_corr(corrs)\n        mean = max(np.mean(np.abs(window1)), np.mean(np.abs(window2)))\n        ts_corr[idx] = corr\n        ts_shift[idx] = shift\n        ts_mean[idx] = mean\n    \n    ts_corr = interpolate_nan(ts_corr)\n    ts_shift = interpolate_nan(ts_shift)\n    ts_mean = interpolate_nan(ts_mean)\n\n    return ts_corr, ts_shift, ts_mean\n\ndef get_corr_anomalies(ts_corr, corr_thres, ts_mean, mean_thres):\n    ts_anomaly = np.logical_and(\n        np.array(ts_corr) < corr_thres, \n        np.array(ts_mean) > mean_thres\n        )\n    return ts_anomaly","metadata":{"execution":{"iopub.status.busy":"2024-02-13T09:53:12.811057Z","iopub.execute_input":"2024-02-13T09:53:12.811396Z","iopub.status.idle":"2024-02-13T09:53:12.833043Z","shell.execute_reply.started":"2024-02-13T09:53:12.811369Z","shell.execute_reply":"2024-02-13T09:53:12.831871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Perform rule-based anomaly detection","metadata":{}},{"cell_type":"code","source":"def apply_rule_based_anomaly_detection(control_area_considered):\n    print(\"Processing control area \" + str(control_area_considered))\n    # Define training and test data sets        \n    df_train = control_area_train_1 if (control_area_considered == 1) else control_area_train_2\n    df_process_list = [df_train]\n    df_process_list_names = ['df_train']\n    df_test = None\n    if (test_data_included):\n        df_test = control_area_test_1 if (control_area_considered == 1) else control_area_test_2\n        for segment_id in df_test['test_data_segment_id'].sort_values().unique():\n            df_process_list.append(df_test.loc[df_test['test_data_segment_id'] == segment_id].copy())\n            df_process_list_names.append('df_test_' + str(segment_id))\n\n    # Anomaly type I: Limit violations for band limited corrected demand\n    print(\"Processing anomaly type I: Limit violations for band limited corrected aFRR demand...\")\n    for df in df_process_list:\n        # Initialize indicator column for limit violations\n        df['flagBandLimitViolation'] = False\n\n        # Detect and mark limit violations\n        df['flagBandLimitViolation'] = (df['BandLimitedCorrectedDemand'] > (df['controlBandPos'] + band_limit_violation_thr)) | \\\n            (df['BandLimitedCorrectedDemand'] < (df['controlBandNeg'] - band_limit_violation_thr))\n\n        # Filter out non participations\n        df['flagBandLimitViolation'] &= (df['participationCMO'] == 1) & (df['participationIN'] == 1)\n\n        # Take only those anomalies into account that last for at least near_const_abs_number_thr time steps\n        df['flagBandLimitViolationFiltered'] = detect_constant_value_ranges(pd.DataFrame(df['flagBandLimitViolation'].astype(int)), 0.1, near_const_abs_number_thr, \\\n                                                ['flagBandLimitViolation'], ['flagBandLimitViolation'], \\\n                                                    fill_first_entry=1, flag_first_const_value=True, diff_first_const_value=1)\n\n    # Anomaly type II: Inconsistencies of demand, correction and corrected demand\n    print(\"Processing anomaly type II: Inconsistencies of aFRR demand, correction and corrected demand...\")\n    for df in df_process_list:\n        # Initialize indicator column for correction inconsistencies\n        df['flagDemandCorrectionInconsistency'] = False\n\n        # Detect and mark inconsistencies\n        df['flagDemandCorrectionInconsistency'] = (df['correctedDemand'] - df['Demand'] - \\\n            df['correction']).abs() > demand_correction_incons_abs_thr\n\n        # Filter out non participations\n        df['flagDemandCorrectionInconsistency'] &= (df['participationCMO'] == 1) & (df['participationIN'] == 1)\n\n        # Take only those anomalies into account that last for at least near_const_abs_number_thr time steps\n        df['flagDemandCorrectionInconsistencyFiltered'] = detect_constant_value_ranges(pd.DataFrame(df['flagDemandCorrectionInconsistency'].astype(int)), 0.1, near_const_abs_number_thr, \\\n                                                    ['flagDemandCorrectionInconsistency'], ['flagDemandCorrectionInconsistency'], \\\n                                                        fill_first_entry=1, flag_first_const_value=True, diff_first_const_value=1)\n\n    # Anomaly type III: Inconsistencies of LFC input + aFRR activation and corrected demand\n    # Initialize indicator column for inconsistencies in LFC input\n    print(\"Processing anomaly type III: Inconsistencies of LFC input + aFRR activation and corrected demand...\")\n    for df in df_process_list:\n        df['flagLFCInputInconsistency'] = False\n    # Calculate absolute deviation threshold\n    aux = (df_train['aFRRactivation'] + df_train['LFCInput'] - df_train['correctedDemand']).abs()\n    LFCInput_incons_thr = aux.quantile(LFCInput_incons_thr_quantile_p)\n    print('LFCInput_incons_thr: ', LFCInput_incons_thr)\n    # Detect and mark inconsistencies\n    for df in df_process_list:\n        aux = (df['aFRRactivation'] + df['LFCInput'] - df['correctedDemand']).abs()\n        df['flagLFCInputInconsistency'] = aux > LFCInput_incons_thr\n\n        # Filter out non participations\n        df['flagLFCInputInconsistency'] &= (df['participationCMO'] == 1) & (df['participationIN'] == 1)\n\n        # Take only those anomalies into account that last for at least near_const_abs_number_thr time steps\n        df['flagLFCInputInconsistencyFiltered'] = detect_constant_value_ranges(pd.DataFrame(df['flagLFCInputInconsistency'].astype(int)), 0.1, near_const_abs_number_thr, \\\n                                                        ['flagLFCInputInconsistency'], ['flagLFCInputInconsistency'], \\\n                                                            fill_first_entry=1, flag_first_const_value=True, diff_first_const_value=1)\n\n    # Anomaly type IV: Extrem values\n    # Define control band width\n    print(\"Processing anomaly type IV: Extrem values...\")\n    for df in df_process_list:\n        df['controlBandWidth'] = df['controlBandPos'] - df['controlBandNeg']\n\n    # Calculate thresholds for negative / positive values as quantiles\n    extreme_values_lower_thrs = df_train.loc[:, select_cols_extreme_values].quantile(extreme_values_lower_thr_quantile_p)\n    extreme_values_upper_thrs = df_train.loc[:, select_cols_extreme_values].quantile(extreme_values_upper_thr_quantile_p)\n    for column in select_cols_extreme_values:\n        if (column in extreme_values_lower_special_thr.index):\n            if (~np.isnan(extreme_values_lower_special_thr[column])):\n                extreme_values_lower_thrs[column] = extreme_values_lower_special_thr[column]\n\n        if (column in extreme_values_upper_special_thr.index):\n            if (~np.isnan(extreme_values_upper_special_thr[column])):\n                extreme_values_upper_thrs[column] = extreme_values_upper_special_thr[column]\n\n    # Initialize indicator columns for extreme values\n    for df in df_process_list:\n        df['flagExtremeValue'] = False\n\n        # Detect and mark extreme values\n        df['flagExtremeValue'] = np.sign((df.loc[:, select_cols_extreme_values] < extreme_values_lower_thrs[select_cols_extreme_values]).sum(axis=1)).astype(bool)\n        df['flagExtremeValue'] |= np.sign((df.loc[:, select_cols_extreme_values] > extreme_values_upper_thrs[select_cols_extreme_values]).sum(axis=1)).astype(bool)\n\n        # Filter out non participations\n        df['flagExtremeValue'] &= (df['participationCMO'] == 1) & (df['participationIN'] == 1)\n\n        # Take only those anomalies into account that last for at least near_const_abs_number_thr time steps\n        df['flagExtremeValueFiltered'] = detect_constant_value_ranges(pd.DataFrame(df['flagExtremeValue'].astype(int)), 0.1, near_const_abs_number_thr, \\\n                                                            ['flagExtremeValue'], ['flagExtremeValue'], \\\n                                                                fill_first_entry=1, flag_first_const_value=True, diff_first_const_value=1)  \n\n    # Anomaly type V: Extreme movements\n    # Initialize indicator columns for extreme jumps\n    print(\"Processing anomaly type V: Extreme movements...\")\n    for df in df_process_list:\n        df['flagExtremeJump'] = False\n\n    # Calculate absolute large movement (\"jump\") thresholds\n    aux = df_train.loc[:, select_cols_extreme_jumps].diff().abs().fillna(0)\n    extreme_jumps_abs_thrs = aux.quantile(extreme_jumps_abs_thr_quantile_p)\n\n    # Detect and mark extreme jumps\n    for df in df_process_list:\n        aux = df.loc[:, select_cols_extreme_jumps].diff().abs().fillna(0)\n        df['flagExtremeJump'] = np.sign((aux > extreme_jumps_abs_thrs[select_cols_extreme_jumps]).sum(axis=1)).astype(bool)\n\n        # Filter out non participations\n        df['flagExtremeJump'] &= (df['participationCMO'] == 1) & (df['participationIN'] == 1)\n\n    # Anomaly type VI: (Near) constant values\n    # Initialize indicator column for (near) constant values\n    print(\"Processing anomaly type VI: (Near) constant values...\")\n    for i, df in enumerate(df_process_list):\n        df['flagNearConstantValue'] = False\n\n        df['flagNearConstantValue'] = detect_constant_value_ranges(df, near_const_abs_thr, near_const_abs_number_thr, \\\n                                                                select_cols_near_const_values, select_cols_near_const_values_elim_0)\n\n        # Filter out non participations\n        df['flagNearConstantValue'] &= (df['participationCMO'] == 1) & (df['participationIN'] == 1)\n\n\n\n    # Anomaly type VII: Inconsistent delays of correction echo\n    # Detect and mark out of phase anomalies between demand correction and correction echo\n    print(\"Processing anomaly type VII: Inconsistent delays of correction echo...\")\n    # Set most frequent shift (to be used if not determined from training data)\n    shift_most_frequent = corr_echo_deviation_shift_most_frequent[control_area_considered - 1]\n    print('shift_most_frequent: ', shift_most_frequent)\n\n    # Determine start index in df_process_list: 0 if training run is to be performed, 1 if not\n    start_index = 1\n    if (corr_echo_incons_delays_train): start_index = 0\n\n    for df_num, df in enumerate(df_process_list[start_index:]):\n        # Initialize indicator column for out of phase shift anomalies\n        df['flagOutOfPhaseCorrection'] = True\n\n        aux = df[['correction', 'correctionEcho']].copy()\n\n        aux_dev = aux.copy()\n        # Generate test shifts and calculate deviation of correction echo and shifted correction\n        for i in range(start_shift_corr, end_shift_corr + 1):\n            # Test shifts\n            aux.loc[:, i] = aux.loc[:, 'correctionEcho'].shift(-i)\n            # Deviations\n            aux_dev.loc[:, i] = (aux.loc[:, 'correction'] - aux.loc[:, i]).abs()\n\n        # Keep only columns with deviations for test shifts\n        aux_dev.drop(columns=['correction', 'correctionEcho'], inplace=True)\n        # Sub threshold deviations encoded by 1, other deviations by 0\n        aux_dev = (aux_dev < corr_echo_deviation_thr).astype(int)\n        # Add extra row for technical reasons (termination of ranges with sub threshold deviations that\n        # extend up to the end of the dataframe by 0 entries)\n        aux_dev = pd.concat([aux_dev, aux_dev[-1:].reindex([aux_dev.index[-1] + 1], fill_value=0)])\n        aux_dev.loc[aux_dev.index[-1:], :] = 0\n\n        # Initialize lists for exemption results\n        shift_list = list()\n        start_index_list = list()\n        end_index_list = list()\n\n        # Perform stepwise exemption of the index range by observation of the expected echo signals\n        # for all possible shifts\n        while True:\n            # Calculate first differences of the 0/1 encoded threshold deviations\n            # Result for each column representing a certain shift:\n            # +1 for the starting row of a consecutive range of matches between the correction and the correction echo\n            # -1 for the row following the consecutive range of matches\n            aux_dev_diff = aux_dev.diff()\n            # Encode first row such that the (unknown) column values for the last row before\n            # the start of the data are assumed not to match \n            aux_dev_diff.loc[aux_dev.index[0], :] = aux_dev.loc[aux_dev.index[0], :]\n\n            max_range_all = None\n\n            # Find the largest matching index range for any of the available shifts\n            # Loop through shifts\n            for i in range(start_shift_corr, end_shift_corr + 1):\n                # Boundaries of matching ranges are characterized by entries of +1 (start) / -1 (end)\n                boundaries = aux_dev_diff.loc[(aux_dev_diff[i] == 1) | (aux_dev_diff[i] == -1)]\n                indexes = boundaries.index\n                # Calculate differences of starting and ending indices to calculate the length of the matching range\n                index_diffs = indexes.to_series().diff()\n                # Proceed if any matching ranges have been found\n                if (len(index_diffs) > 0):\n                    # Calculate maximum length of the matching ranges for the current shift\n                    # Take only into account differences ending with a -1 entry in the deviation dataframe\n                    # (end of each matching range)\n                    max_range = index_diffs[boundaries[i] == -1].max()\n                    if ((max_range_all is None) or (max_range > max_range_all)):\n                        # Save data for current maximum matching range if larger than the ranges found\n                        # for the previous shifts in the loop\n                        shift_max_range_all = i\n                        max_range_all = max_range\n                        max_range_arg_all = index_diffs[boundaries[i] == -1].argmax()\n                        end_index_max_range_arg_all = indexes[boundaries[i] == -1][max_range_arg_all]\n                        start_index_max_range_arg_all = indexes[indexes.get_loc(end_index_max_range_arg_all) - 1]\n                        end_index_max_range_arg_all = end_index_max_range_arg_all - 1\n\n            # Check if any further matching ranges have been found\n            if(max_range_all is not None):\n\n                # Apply exemption for the largest macthing range found, i.e. put a code of 0 in the deviation dataframe\n                # to mask that matching range from any further processing \n                aux_dev.loc[start_index_max_range_arg_all:end_index_max_range_arg_all, :] = 0\n\n                # Update exemption lists\n                shift_list.append(shift_max_range_all)\n                start_index_list.append(start_index_max_range_arg_all)\n                end_index_list.append(end_index_max_range_arg_all)\n\n            else:\n                # Break if no further matching ranges have been found\n                break\n\n        aux = df[['correction', 'correctionEcho']].copy()\n        aux['correctionEchoShift'] = np.nan\n        aux['correctionTest'] = np.nan\n        aux['referencedIndex'] = -1\n\n        for i in range(len(shift_list)):\n\n            # Calculate the echo shift in time steps for each matching range found during the exemption processing\n            aux.loc[start_index_list[i]:end_index_list[i], 'correctionEchoShift'] = shift_list[i]\n            # Calculate the deviation of correction and correctionEcho shifted by the echo shift determined before for validation\n            aux.loc[start_index_list[i]:end_index_list[i], 'correctionTest'] = aux.loc[start_index_list[i]:end_index_list[i], 'correction'] - \\\n                aux.shift(-shift_list[i]).loc[start_index_list[i]:end_index_list[i], 'correctionEcho']\n\n            # Set start and end location of the shifted index\n            start_loc = (aux.index.get_loc(start_index_list[i]) + shift_list[i])\n            end_loc = (aux.index.get_loc(end_index_list[i]) + shift_list[i])\n            # Calculate the index referred to by the current index + shift\n            aux.loc[start_index_list[i]:end_index_list[i], 'referencedIndex'] = aux.index[start_loc:(end_loc + 1)]\n\n        aux['idx'] = aux.index\n        aux = aux.sort_values(['referencedIndex', 'idx'])\n        aux['idxDiff'] = aux['referencedIndex'].diff()\n        aux.loc[aux['idxDiff'] == 0, 'correctionEchoShift'] = np.nan\n        aux.sort_index(inplace=True)\n\n        print('loop, dataframe, sum correctionTest, correctionEchoShift unique: ', df_num, df_process_list_names[df_num],  \\\n              aux['correctionTest'].sum(), aux['correctionEchoShift'].unique())\n\n        if (corr_echo_incons_delays_train and (df_num == 0)):\n            # Determine the shift with the largest number of sub threshold deviations for the training data\n            shift_most_frequent = int(aux.loc[:, 'correctionEchoShift'].value_counts().idxmax())\n            print(aux.loc[:, 'correctionEchoShift'].value_counts(dropna=False))\n            print('shift_most_frequent: ', shift_most_frequent)\n\n        df.loc[df.index, 'correctionEchoShift'] = aux.loc[df.index, 'correctionEchoShift']\n\n        # No out of phase anomaly is detected for sub threshold deviations of the current shift from the most frequent shift\n        df.loc[(df['correctionEchoShift'] - shift_most_frequent).abs() <= corr_echo_deviation_shift_most_frequent_thr, 'flagOutOfPhaseCorrection'] = False\n        # Shift deviations for the correction echo smaller than 1 should not occur since the echo signal should at least take one timestep\n        df.loc[df['correctionEchoShift'] < 1, 'flagOutOfPhaseCorrection'] = True\n        # An anomaly is detected if a correction signal could not identified in the correction echo\n        df.loc[df['correctionEchoShift'].isna(), 'flagOutOfPhaseCorrection'] = True\n\n        # No out of phase anomaly is detected for the last time steps determined by the most frequent shift\n        # as the corresponding correction values are not present in the data (correctionEcho)\n        df.loc[df.index[-shift_most_frequent:], 'flagOutOfPhaseCorrection'] = False\n\n        # Filter out non participations\n        df['flagOutOfPhaseCorrection'] &= (df['participationCMO'] == 1) & (df['participationIN'] == 1)\n\n    # Free memory\n    del aux\n    del aux_dev\n    # Trigger garbage collection\n    gc.collect()\n\n\n    # Anomaly type VIII: Asynchronous activation of aFRR\n    # Calculate simulated aFRR activation\n    print(\"Processing anomaly type VIII: Asynchronous activation of aFRR...\")\n    # Control area 1: The optimal result of 81 time steps for the training set of control area 1 is only slightly above the expected 75 time steps according to regulatory requirements\n    # Control area 2: The optimal 16 time steps indicate that aFRRactivation usually follows aFRR request more tightly than in control area 1\n    for df in df_process_list:\n        window_size = window_size_aFRRactivationSim[control_area_considered - 1]\n        response_time = 0\n        df['aFRRactivationSim'] = df['aFRRrequest'].rolling(window_size - response_time, min_periods=1).mean().shift(response_time)\n\n    # Find time ranges with large absolute differences between actual aFRR activation and simulated aFRR activation / aFRR request\n\n    ae_abs_thr = None\n    ae_jump_abs_thr = None\n\n    for df_num, df in enumerate(df_process_list):\n        print()\n        print('Dataframe: ', df_process_list_names[df_num])\n\n        print('Number of rows: ', len(df))\n\n        # Calculate distance for actual activation\n        dist_req = df['aFRRrequest'] - df['aFRRactivation']\n\n        # Calculate distance for simulated activation\n        dist_sim = df['aFRRactivationSim'] - df['aFRRactivation']\n        dist_sim.loc[dist_sim.isna()] = np.inf\n\n        # Calculate minimum absolute distances between actual activation and simulated activation / actual request \n        ae = pd.Series(np.where(dist_sim.abs() < dist_req.abs(), dist_sim, dist_req))\n        ae.index = df.index\n\n        # Calculate running differences for the distance\n        ae_jump_abs = ae.diff().fillna(0).abs()\n\n        if (df_num == 0):\n            # Calculate threshold value for the distance movement as a quantile (only for training data)\n            ae_jump_abs_thr = ae_jump_abs.quantile(ae_jump_abs_thr_quantile_p)\n            print('ae_jump_thr: ', ae_jump_abs_thr)\n            # Calculate threshold as quantile for the absolute distance (only for training data)\n            ae_abs_thr = ae.abs().quantile(ae_thr_quantile_p)\n            print('ae_thr: ', ae_abs_thr)\n\n        # Peform further calculations for test data\n        # Run with training data is included depending on parameter setting\n        if (ae_include_train or (df_num > 0)):\n\n            # Set index values for absolute distances above threshold\n            # No anomaly if the distance of aFRRactivation from aFRRrequest and aFRRactivationSim has different signs, i.e. aFRRactivation is\n            # located between aFRRrequest and aFRRactivationSim\n            ae_above_thr_idx = pd.Series(df[(ae.abs() > ae_abs_thr) & (np.sign(dist_req) == np.sign(dist_sim))].index.values)\n\n            # Detect lower and upper bounds for ranges with absolute distances above threshold\n            ae_above_thr_bounds_lower_idx = (ae_above_thr_idx - ae_above_thr_idx.shift(1, fill_value=np.inf)).abs() > 1\n            ae_above_thr_bounds_upper_idx = (ae_above_thr_idx - ae_above_thr_idx.shift(-1, fill_value=np.inf)).abs() > 1\n\n            # Store lower and upper bounds in dataframe\n            ae_above_thr_bounds_idx = pd.DataFrame(data={ \\\n            'lower_bound_thr': ae_above_thr_idx[ae_above_thr_bounds_lower_idx].values, \\\n            'upper_bound_thr': ae_above_thr_idx[ae_above_thr_bounds_upper_idx].values})\n\n            # Detect index locations where the actual activation comes close either to the simulated activation or (and) the request (\"near locations\")\n\n            # Detect index locations where the minimum distance is below threshold\n            ae_near_idx = pd.DataFrame(data={'idx': ae[ae.abs() < ae_near_abs_thr].index})\n\n            # Detect which near locations qualify for the start or end of an anomaly\n            # For this, the two movements of the distance adjacent to the left (start of anomaly) or right (end of anomaly)\n            # have to show moderate values (i.e. no \"jumps\")\n\n            # Detect index values with distance movements below threshold\n            ae_jump_abs_below_thr = ae_jump_abs < ae_jump_abs_thr\n\n            # Qualification for start / end of anomaly requires two adjacent movements of activation error below threshold to the left / right\n            anomaly_qualified_start = ae_jump_abs_below_thr & ae_jump_abs_below_thr.shift(1, fill_value=True)\n            anomaly_qualified_end = ae_jump_abs_below_thr.shift(-1, fill_value=True) & ae_jump_abs_below_thr.shift(-2, fill_value=True)\n\n            # Store information regarding qualification with the near indexes\n            ae_near_idx.loc[:, 'anomaly_qualified_start'] = anomaly_qualified_start[ae_near_idx.idx].values\n            ae_near_idx.loc[:, 'anomaly_qualified_end'] = anomaly_qualified_end[ae_near_idx.idx].values\n\n            # Set start and end values of index ranges that represent anomalies (\"anomaly windows\")\n\n            # Init new columns to ensure integer values\n            ae_above_thr_bounds_idx['start_anom'] = -1\n            ae_above_thr_bounds_idx['end_anom'] = -1\n\n            upper_bound_index = max(ae_jump_abs.index)\n\n            # For each anomaly window, set start / end values using the nearest qualified sub-threshold distance points to the left / right\n            ae_above_thr_bounds_idx.loc[:, 'start_anom'] = [min(np.max(ae_near_idx.loc[(ae_near_idx['idx'] <= ae_above_thr_bounds_idx.loc[i, 'lower_bound_thr']) & \\\n                ae_near_idx['anomaly_qualified_start'], 'idx']) + 1, upper_bound_index) for i in ae_above_thr_bounds_idx.index]\n\n            ae_above_thr_bounds_idx.loc[ae_above_thr_bounds_idx['start_anom'].isna(), 'start_anom'] = \\\n                ae_above_thr_bounds_idx.loc[ae_above_thr_bounds_idx['start_anom'].isna(), 'lower_bound_thr']\n\n            ae_above_thr_bounds_idx.loc[:, 'end_anom'] = [max(np.min(ae_near_idx.loc[(ae_near_idx['idx'] >= ae_above_thr_bounds_idx.loc[i, 'upper_bound_thr']) & \\\n                ae_near_idx['anomaly_qualified_end'], 'idx']) - 1, 0) for i in ae_above_thr_bounds_idx.index]\n\n            ae_above_thr_bounds_idx.loc[ae_above_thr_bounds_idx['end_anom'].isna(), 'end_anom'] = \\\n                ae_above_thr_bounds_idx.loc[ae_above_thr_bounds_idx['end_anom'].isna(), 'upper_bound_thr']\n\n            # Provide indicator column for detected anomalies\n            df['flagAsynchronousaFRRactivation'] = False\n\n            for i in ae_above_thr_bounds_idx.index:\n                df.loc[ae_above_thr_bounds_idx.loc[i, 'start_anom']:ae_above_thr_bounds_idx.loc[i, 'end_anom'], \\\n                'flagAsynchronousaFRRactivation'] = True\n\n            # Filter out non participations\n            df['flagAsynchronousaFRRactivation'] &= (df['participationCMO'] == 1) & (df['participationIN'] == 1)\n\n            # Check results\n\n            # Set margin of index values to the left and right of anomaly windows\n            margin = 1\n            aux = ae_above_thr_bounds_idx[['start_anom', 'end_anom']].drop_duplicates(keep='first')\n\n\n    # Anomaly type IX: Shifted correlation of aFRRactivation and aFRRrequest\n    print(\"Processing anomaly type IX: Shifted correlation of aFRRactivation and aFRRrequest...\")\n    for df in df_process_list[1:]:\n        timeseries1 = df[\"aFRRrequest\"]\n        timeseries2 = df[\"aFRRactivation\"]\n        ts_corr, ts_shift, ts_mean = get_shifed_corr_series(timeseries1, timeseries2, window_size_corr, min_shift_corr, max_shift_corr)\n        ts_anomaly = get_corr_anomalies(ts_corr, corr_thr, ts_mean, corr_noise_thr)\n        df[\"flagShiftedCorrelationRequestActivation\"] = ts_anomaly\n        \n        \n    # Put together results for test data segments\n    if (test_data_included):\n        df_test = pd.concat(df_process_list[1:])\n        \n    for df in [df_train, df_test] if test_data_included else [df_train]:\n        aux = df[[column for column in df.columns if column.startswith('flag')]].sum(axis=1)\n        df['flagAnomalyRuleBased'] = aux > 0\n        \n    aux = df_test\n    # aux = df_train\n    \n    # Print statistics\n    for column in [column for column in aux.columns if (column.startswith('flag') | column.startswith('1flag'))]:\n        perc = len(aux[aux[column] == True]) / len(aux)\n        print('flag, perc.: ' + column +', %8.4f' % (perc))\n    \n    df_train.to_csv(f\"processed_area{control_area_considered}_train.csv\", index=True)\n    df_test.to_csv(f\"processed_area{control_area_considered}_test.csv\", index=True)","metadata":{"execution":{"iopub.status.busy":"2024-02-13T09:54:54.651474Z","iopub.status.idle":"2024-02-13T09:54:54.651994Z","shell.execute_reply.started":"2024-02-13T09:54:54.651739Z","shell.execute_reply":"2024-02-13T09:54:54.651760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Split input data according to control areas \n# Training set\ncontrol_area_train_1 = raw_train_data[raw_train_data['controlArea']==1].copy()\ncontrol_area_train_2 = raw_train_data[raw_train_data['controlArea']==2].copy()\n\"\"\"\n'Datum_Uhrzeit_CET', 'Demand', 'correction', 'correctedDemand', 'FRCE',\n       'controlBandPos', 'controlBandNeg', 'LFCInput', 'aFRRactivation',\n       'aFRRrequest', 'participationCMO', 'participationIN', 'correctionEcho',\n       'BandLimitedCorrectedDemand', 'controlArea'\n\"\"\"\n\n# Test set\ncontrol_area_test_1 = raw_test_data[raw_test_data['controlArea']==1].copy()\ncontrol_area_test_2 = raw_test_data[raw_test_data['controlArea']==2].copy()\n\"\"\"\n'Datum_Uhrzeit_CET', 'Demand', 'correction', 'correctedDemand',\n       'FRCE', 'controlBandPos', 'controlBandNeg', 'LFCInput',\n       'aFRRactivation', 'aFRRrequest', 'participationCMO', 'participationIN',\n       'correctionEcho', 'BandLimitedCorrectedDemand', 'controlArea',\n       'test_data_segment_id'\n\"\"\"\n\n# Free memory\ndel raw_train_data\n# Trigger garbage collection\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-02-13T09:54:54.653266Z","iopub.status.idle":"2024-02-13T09:54:54.653725Z","shell.execute_reply.started":"2024-02-13T09:54:54.653495Z","shell.execute_reply":"2024-02-13T09:54:54.653515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_data_included = True\n\nfor control_area_considered in [1, 2]:\n    apply_rule_based_anomaly_detection(control_area_considered)","metadata":{"execution":{"iopub.status.busy":"2024-02-13T09:54:55.551175Z","iopub.execute_input":"2024-02-13T09:54:55.551544Z","iopub.status.idle":"2024-02-13T09:58:26.597678Z","shell.execute_reply.started":"2024-02-13T09:54:55.551511Z","shell.execute_reply":"2024-02-13T09:58:26.595994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Merge ControlAreas","metadata":{}},{"cell_type":"code","source":"rb_train_area1 = pd.read_csv(\"/kaggle/working/processed_area1_train.csv\", parse_dates=['Datum_Uhrzeit_CET']).sort_index()\nrb_train_area2 = pd.read_csv(\"/kaggle/working/processed_area2_train.csv\", parse_dates=['Datum_Uhrzeit_CET']).sort_index()\n\nrb_test_area1 = pd.read_csv(\"/kaggle/working/processed_area1_test.csv\", parse_dates=['Datum_Uhrzeit_CET']).sort_index()\nrb_test_area2 = pd.read_csv(\"/kaggle/working/processed_area2_test.csv\", parse_dates=['Datum_Uhrzeit_CET']).sort_index()\nrb_test = pd.concat([rb_test_area1, rb_test_area2])","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:30:29.002081Z","iopub.execute_input":"2024-02-13T10:30:29.002486Z","iopub.status.idle":"2024-02-13T10:31:01.712729Z","shell.execute_reply.started":"2024-02-13T10:30:29.002452Z","shell.execute_reply":"2024-02-13T10:31:01.711625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Machine Learning approach","metadata":{}},{"cell_type":"markdown","source":"## Preprocessing\n### Filtering of columns\nNot all columns are equally relevant for a good prediction. We predefine them to later loop over the relevant ones. ","metadata":{}},{"cell_type":"code","source":"relevant_cols = [\n    'Datum_Uhrzeit_CET',\n    'Demand',\n    'correction',\n    'correctedDemand',\n    'FRCE',\n    'LFCInput',\n    'aFRRactivation',\n    'aFRRrequest',\n    'correctionEcho',\n    'BandLimitedCorrectedDemand',\n    'controlArea'\n    ]\nvalue_cols = [col for col in relevant_cols if col not in [\"Datum_Uhrzeit_CET\",\"controlArea\"]]\nhelp_cols = ['controlBandPos', 'controlBandNeg']\nX_train = X_train[relevant_cols + help_cols]\nX_test_filter = X_test[[\"participationIN\", \"participationCMO\"]]\nX_test = X_test[relevant_cols+help_cols + [\"test_data_segment_id\"]]","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:08:18.687124Z","iopub.execute_input":"2024-02-13T10:08:18.687407Z","iopub.status.idle":"2024-02-13T10:08:18.919694Z","shell.execute_reply.started":"2024-02-13T10:08:18.687383Z","shell.execute_reply":"2024-02-13T10:08:18.918849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Collection of feature engineering functions","metadata":{}},{"cell_type":"code","source":"def prepare_shifts(df, value_cols, shifts):\n    # shift data in terms of index steps\n    for shift in shifts:\n        df[[col + f\"_shift_{shift}\" for col in value_cols]] = df[value_cols].shift(shift, fill_value=0)\n    return df\n\ndef prepare_diffs(df, value_cols, diffs):\n    # calculate diff of shifted and original data\n    for diff in diffs:\n        df[[col + f\"_diff_{diff}\" for col in value_cols]] = df[value_cols].shift(diff, fill_value=0) - df[value_cols]\n    return df\n\ndef prepare_rolling_window(df, value_cols, periods):\n    # calculate rolling mean\n    for period in periods:\n        df[[col + f\"_rolling_mean_{period}\" for col in value_cols]] = df[value_cols].rolling(period).mean().fillna(method='ffill').fillna(method='bfill')\n    return df\n\ndef prepare_rolling_window_std(df, value_cols, periods):\n    # calculate rolling std deviation\n    for period in periods:\n        df[[col + f\"_rolling_std_{period}\" for col in value_cols]] = df[value_cols].rolling(period).std().fillna(method='ffill').fillna(method='bfill')\n    return df\n\ndef prepare_control_bands(df):\n    # define control bands\n    diff_band_pos = df[\"controlBandPos\"] - df[\"aFRRactivation\"]\n    diff_band_neg = df[\"aFRRactivation\"] - df[\"controlBandNeg\"]\n    df[\"diff_band\"] = np.minimum(diff_band_pos, diff_band_neg)\n    df[\"controlBandWindow\"] = df[\"controlBandPos\"] - df[\"controlBandNeg\"]\n    df = df.drop(columns=[\"controlBandPos\",\"controlBandNeg\"], axis=1)\n    return df\n\ndef prepare_time(df):\n    # turned out to have no relevance, especially, as training data are shifted in time\n    df['Datum_Uhrzeit_CET'] = pd.to_datetime(df['Datum_Uhrzeit_CET'])\n\n    # Create new features from the timestamp\n    # df['year'] = df['Datum_Uhrzeit_CET'].dt.year\n    # df['month'] = df['Datum_Uhrzeit_CET'].dt.month\n    df['day'] = df['Datum_Uhrzeit_CET'].dt.day\n    df['hour'] = df['Datum_Uhrzeit_CET'].dt.hour\n    df['minute'] = df['Datum_Uhrzeit_CET'].dt.minute\n    df['second'] = df['Datum_Uhrzeit_CET'].dt.second\n    df['day_of_week'] = df['Datum_Uhrzeit_CET'].dt.dayofweek  # Monday=0, Sunday=6.\n    df['is_weekend'] = df['day_of_week'].isin([5, 6]).astype(int)  # Create a boolean column: True for weekend.\n\n    # df.drop(\"Datum_Uhrzeit_CET\", axis=1)\n    return df\n\ndef filter_unimportant_data(df):\n    # filter out data without participation flag\n    df = df.loc[df[\"participationCMO\"]==1 &\tdf[\"participationIN\"]==1]\n    return df\n\ndef prepare_correlation(df,corr_value_cols, window):\n    # definition of rolling correlation\n    for (var1,var2) in corr_value_cols:\n        df[f\"{var1}_{var2}_corr_{window}\"]=df[var1].rolling(window=window).corr(df[var2]).fillna(0)\n    return df\n","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:08:22.297132Z","iopub.execute_input":"2024-02-13T10:08:22.297501Z","iopub.status.idle":"2024-02-13T10:08:22.313676Z","shell.execute_reply.started":"2024-02-13T10:08:22.297474Z","shell.execute_reply":"2024-02-13T10:08:22.312554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Preprocessing main script","metadata":{}},{"cell_type":"code","source":"def preprocessing(df, test = False, scaler = None):\n\n    # general preprocessing for all useful variables:\n    # The values for rolling windows and diffs were mainly found by testing and technical explanations\n    df = prepare_rolling_window(df, value_cols, [10, 40])\n    smooth_value_cols = [col + \"_rolling_mean_10\" for col in value_cols]\n    df = prepare_diffs(df, smooth_value_cols, [15, 40])\n    df = prepare_diffs(df, value_cols, [1])\n    \n    # shifts sometimes seemed to lead to errors, so they are not used:\n    # df = prepare_shifts(df, smooth_value_cols, [15, 40])\n\n    # Calculate noise and first differences as a derivative approximation\n    for col in value_cols:\n        df[col + \"_noise\"] = df[col] - df[col + \"_rolling_mean_10\"]\n        df[col + \"_short_diff_1\"] = df[col] - df[col + \"_diff_1\"]\n\n    # Try to premodel the connection between activation, request and other interdepending variables\n    # We model the correlation and the differences over time\n    df = prepare_shifts(df,[\"aFRRrequest_rolling_mean_10\"],[40])\n    df = prepare_shifts(df,[\"correctedDemand_rolling_mean_10\"],[75])\n    df = prepare_correlation(df,[(\"aFRRrequest_rolling_mean_10_shift_40\",\"aFRRactivation_rolling_mean_10\"),\n                                  (\"aFRRrequest_rolling_mean_10\",\"correctedDemand_rolling_mean_10\"),\n                                  (\"aFRRactivation_rolling_mean_10\",\"correctedDemand_rolling_mean_10\"),\n                                  (\"FRCE_rolling_mean_10\",\"correctedDemand_rolling_mean_10\")],\n                             100)\n\n    df[\"aFRRrequest_rolling_mean_10_shift_40_aFRRactivation_rolling_mean_10_diff\"] = df[\"aFRRrequest_rolling_mean_10_shift_40\"]-df[\"aFRRactivation_rolling_mean_10\"]\n    df[\"aFRRactivation_rolling_mean_10_correctedDemand_rolling_mean_10_shift_75\"] = df[\"aFRRactivation_rolling_mean_10\"] - df[\"correctedDemand_rolling_mean_10_shift_75\"]\n\n    # corrected demand manually\n    df[\"correctedDemand_calc\"] = df[\"Demand\"] + df[\"correction\"]\n    df[\"Demand_Delta\"] = df[\"correctedDemand\"] - df[\"correctedDemand_calc\"]\n    df = prepare_rolling_window(df,[\"Demand_Delta\"], [15])\n\n    # Calculate the distance to the control band limits\n    df = prepare_control_bands(df)\n\n    # Select columns that should be scaled (all that are actual relevant for a prediction)\n    scaler_cols = [col for col in df.columns if col not in [\"Datum_Uhrzeit_CET\",\"controlArea\",\"test_data_segment_id\"]]\n\n    # replace and mark infinity values in the correlation features (e.g. by constant intervals)\n    correlation_cols = [col for col in df.columns if \"_corr_\" in col]\n    df['inf_present'] = df[correlation_cols].apply(lambda row: int(np.isinf(row).any()), axis=1)\n    df = df.replace(to_replace=[np.inf, -np.inf], value=0)\n\n    # Fit and scale the whole train set, if input is the testset then don't fit, just apply the fitted scaler\n    # Use seperate scalers for each controlArea \n    if not test:\n        scaler1 = StandardScaler()\n        scaler2 = StandardScaler()\n        if sum(df[\"controlArea\"] == 1)>0:\n            df.loc[df[\"controlArea\"] == 1, scaler_cols] = scaler1.fit_transform(df.loc[df[\"controlArea\"] == 1, scaler_cols])\n        if sum(df[\"controlArea\"] == 2)>0:\n            df.loc[df[\"controlArea\"] == 2, scaler_cols] = scaler2.fit_transform(df.loc[df[\"controlArea\"] == 2, scaler_cols])\n        return df, (scaler1, scaler2)\n\n    else:\n        if sum(df[\"controlArea\"] == 1)>0:\n            df.loc[df[\"controlArea\"] == 1, scaler_cols] = scaler[0].transform(df.loc[df[\"controlArea\"] == 1, scaler_cols])\n        if sum(df[\"controlArea\"] == 2)>0:\n            df.loc[df[\"controlArea\"] == 2, scaler_cols] = scaler[1].transform(df.loc[df[\"controlArea\"] == 2, scaler_cols])\n        return df","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:08:23.329348Z","iopub.execute_input":"2024-02-13T10:08:23.330016Z","iopub.status.idle":"2024-02-13T10:08:23.346216Z","shell.execute_reply.started":"2024-02-13T10:08:23.329985Z","shell.execute_reply":"2024-02-13T10:08:23.345191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def prepare_testset(X_test, scaler1 = None, scaler2 = None):\n    df_test = pd.DataFrame()\n    for segment in X_test.test_data_segment_id.unique():\n        df_test = pd.concat([df_test, preprocessing(X_test.loc[X_test[\"test_data_segment_id\"]==segment], test=True, scaler=(scaler1,scaler2))], ignore_index=True)\n    return df_test","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:08:23.780899Z","iopub.execute_input":"2024-02-13T10:08:23.781793Z","iopub.status.idle":"2024-02-13T10:08:23.787088Z","shell.execute_reply.started":"2024-02-13T10:08:23.781760Z","shell.execute_reply":"2024-02-13T10:08:23.786135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train_1 = X_train.loc[X_train[\"controlArea\"] == 1]\nX_train_2 = X_train.loc[X_train[\"controlArea\"] == 2]\nX_train_1 = X_train_1[::4]\nX_train_2 = X_train_2[::4]\ndf_1, (scaler1,_) = preprocessing(X_train_1)\ndf_2, (_,scaler2) = preprocessing(X_train_2)","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:08:24.392820Z","iopub.execute_input":"2024-02-13T10:08:24.393165Z","iopub.status.idle":"2024-02-13T10:11:06.607197Z","shell.execute_reply.started":"2024-02-13T10:08:24.393140Z","shell.execute_reply":"2024-02-13T10:11:06.606152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test = prepare_testset(X_test, scaler1, scaler2)","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:38:34.420232Z","iopub.execute_input":"2024-02-13T10:38:34.421181Z","iopub.status.idle":"2024-02-13T10:38:38.884808Z","shell.execute_reply.started":"2024-02-13T10:38:34.421141Z","shell.execute_reply":"2024-02-13T10:38:38.883774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Use rule-based anomalies as features for the machine learning approach","metadata":{}},{"cell_type":"code","source":"rb_columns = [c for c in rb_train_area1.columns if \"flag\" in c]","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:12:35.333264Z","iopub.execute_input":"2024-02-13T10:12:35.333675Z","iopub.status.idle":"2024-02-13T10:12:35.338664Z","shell.execute_reply.started":"2024-02-13T10:12:35.333646Z","shell.execute_reply":"2024-02-13T10:12:35.337616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rb_train_area1[rb_columns] = rb_train_area1[rb_columns].astype(int)\nrb_train_area2[rb_columns] = rb_train_area2[rb_columns].astype(int)\n\ndf_1 = pd.merge(df_1, rb_train_area1[[\"Datum_Uhrzeit_CET\"] + rb_columns], how=\"left\", on=\"Datum_Uhrzeit_CET\")\ndf_2 = pd.merge(df_2, rb_train_area2[[\"Datum_Uhrzeit_CET\"]+ rb_columns], how=\"left\", on=\"Datum_Uhrzeit_CET\")","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:31:01.714745Z","iopub.execute_input":"2024-02-13T10:31:01.715525Z","iopub.status.idle":"2024-02-13T10:31:04.866349Z","shell.execute_reply.started":"2024-02-13T10:31:01.715487Z","shell.execute_reply":"2024-02-13T10:31:04.865310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rb_test[rb_columns] = rb_test[rb_columns].astype(int)\ndf_test = pd.merge(df_test, rb_test[[\"Datum_Uhrzeit_CET\", \"test_data_segment_id\"]+ rb_columns], how=\"left\", on=[\"Datum_Uhrzeit_CET\", \"test_data_segment_id\"])","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:38:46.641606Z","iopub.execute_input":"2024-02-13T10:38:46.641964Z","iopub.status.idle":"2024-02-13T10:38:46.701018Z","shell.execute_reply.started":"2024-02-13T10:38:46.641935Z","shell.execute_reply":"2024-02-13T10:38:46.700044Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_1.to_csv(f\"final_area1_train.csv\", index=False)\ndf_2.to_csv(f\"final_area2_train.csv\", index=False)\ndf_test.to_csv(f\"final_test.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2024-02-13T10:39:14.024773Z","iopub.execute_input":"2024-02-13T10:39:14.025494Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Time-series LSTM-based Variational Autoencoder Model\n\nA Time-series Long Short-Term Memory (LSTM) based Variational Autoencoder (VAE) model utilizes the power of both LSTMs and VAEs, and is specially designed for time-series data analysis. Here's a simplified explanation of how it works:\n\n1. Understanding the Time-series data & LSTM: Time-series data is essentially a collection of observations made sequentially over time. A LSTM is a type of recurrent neural network (RNN) that is able to learn long-term dependencies, which is particularly important when working with time-series data, as the understanding of a current point in time might depend on a few points (or a ‘sequence’) in the past.\n1. Understanding the Variational Autoencoder (VAE): An autoencoder is a type of artificial neural network used for learning efficient codings of input data. It is designed to reconstruct input data by first encoding the input into a lower dimensional latent-space representation, and then decoding this representation back to the original input. A Variational Autoencoder (VAE), unlike a traditional autoencoder, makes the assumption that the latent space should model a probability distribution (Gaussian). Therefore the VAE learns to generate parameters (mean and variance) of a distribution from which the latent representation of the input data is sampled. This sampled representation is then decoded back to the original input space. The main goal of a VAE is learning to reconstruct the input data.\n1. Combining these concepts into a Time-series LSTM-based VAE model: In a time-series LSTM-based VAE model, the encoder part of the VAE is designed using LSTM layers to better learn the interdependencies of the time-series data. The decoded outputs are also time-series values which are generated from the latent space.\nThe workflow of such a model can be summarized in four steps:\n    - Encoding: The LSTM-based encoder processes the time-series data of a specific window size and generates a dense representation. This model is designed to encode temporal dependencies of time-series patterns.\n    - Gaussian Sampling: The dense representation is then used to parameterize a Gaussian probability distribution, i.e., by learning the mean and the standard deviation.\n    - Decoding: The decoder, also usually an LSTM, takes random samples from the Gaussian distribution generated by the encoder and attempts to reconstruct the original time-series data.\n    - Training: The difference between the original and the reconstructed data forms the reconstruction loss.\nThis approach is commonly used in anomaly detection in time-series data, as anomalies are usually transcribed as data points that the model finds difficult to reconstruct. In addition we add a KL-loss component, which forces the distribution generated by the encoder to be similar to the prior probability of the input vector, which is assumed to be normal, and therefore pushes the latent feature space towards normality.","metadata":{}},{"cell_type":"markdown","source":"### Defintion of utils-functions","metadata":{}},{"cell_type":"code","source":"def preprocess_for_autoencoder(element):\n    # Autoencoder input and output are the same\n    return (element, element)\n\ndef reconstruct_data_segment(df, window_size, model):\n\n    # initialization of dataframe and array\n    df_reconstructed = pd.DataFrame(columns=df.columns)\n    \n\n    for row in range(df.shape[0] - window_size +1):\n        df_sample = df[row : row + window_size]\n        indeces = df_sample.index\n\n        array_truth = df_sample.values[np.newaxis,:,:]\n\n        # let model predict the sample\n        array_pred = model.predict(array_truth, verbose=0)\n\n        # save as array\n        #array_reconstructed[row] = np.squeeze(array_pred)\n\n        # save as dataframe\n        df_sample_reconstructed = pd.DataFrame(data = np.squeeze(array_pred), index =indeces , columns=df_sample.columns)\n\n        # concat dataframe\n        df_reconstructed = pd.concat([df_reconstructed, df_sample_reconstructed])\n\n    # mean prediction of each time step\n    df_reconstructed_mean = df_reconstructed.groupby(level=0).mean()\n\n    return df_reconstructed_mean, df_reconstructed\n\ndef recontruction_to_anomalies(df, df_reconstructed):\n\n    mse = ((df_reconstructed - df )**2 ).sum(axis = 1)\n    return pd.DataFrame(mse, columns = ['anomaly'])\n\ndef create_dataset(self, path, group_ids, slice_length, stride):\n    self.num_features = len(pd.read_csv(path,nrows=1).columns) - 1\n    column_defaults =  [tf.float32]*self.num_features\n    print(self.num_features)\n\n    dataset = tf.data.experimental.CsvDataset(\n        filenames=path,\n        record_defaults=column_defaults,\n        header=True,\n        exclude_cols=[0]\n    )\n    dataset = dataset.map(lambda *cols: tf.stack(cols, axis=0))\n\n    # List to hold datasets for each group\n    group_datasets = []\n\n    # For each group, slice the data corresponding to the group and\n    # create a windowed dataset. Then add it to the list.\n    for group in np.unique(group_ids):\n        indices = [i for i, x in enumerate(group_ids) if x == group]\n        group_dataset = tf.data.Dataset.from_tensor_slices(list(dataset.as_numpy_iterator())[min(indices):max(indices)+1])\n        group_dataset = group_dataset.window(size=slice_length, shift=stride, drop_remainder=True)\n        group_dataset = group_dataset.flat_map(lambda *window: tf.data.Dataset.zip(window).batch(slice_length))\n        group_datasets.append(group_dataset)\n\n    # Concatenate datasets of all groups\n    dataset = group_datasets[0]\n    for gds in group_datasets[1:]:\n        dataset = dataset.concatenate(gds)\n\n    return dataset","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Model definition","metadata":{}},{"cell_type":"code","source":"class TimeSeries_LSTM_Variational_Autoencoder():\n    def __init__(self,\n                 model=None,\n                #  window_size=16,\n                 threshold=None,\n                 q=0.8):\n        \n        # self.controlArea = controlArea\n        self.model = model # tensorflow model\n        # self.slice_length = window_size # input sequence length of the model\n        self.threshold = threshold # threshold to predict anomalies\n        self.threshold_1 = None # threshold to predict anomalies on area 1\n        self.threshold_2 = None # threshold to predict anomalies on area 2\n        self.q = q # quantil values, that is used to generate a threshold when no threshold is given, e.g. q = 0.7 -> we identify an anomalie if the mse of one sample is greater than 70% of the mses of the whole control area\n\n        self.all_anomalies = None\n        self.area_1_anomalies = None\n        self.area_2_anomalies = None\n\n        self.all_reconstructed = None\n        self.area_1_reconstructed = None\n        self.area_2_reconstructed = None\n\n    def configure_for_training(self,\n                               path_area1,\n                               path_area2,\n                               skipping=1,\n                               slice_length=32,\n                               batch_size=64*4):\n        \n        self.slice_length = slice_length\n        self.batch_size = batch_size\n        self.path_area1 = path_area1\n        self.path_area2 = path_area2\n        self.skipping = skipping\n\n        self.num_features = len(pd.read_csv(path_area1,nrows=1).columns) - 1\n        column_defaults =  [tf.float32]*self.num_features\n\n        # steps per epoch, batch size and slice length should not be larger than number of datapoints to not overfit\n        self.steps_per_epoch = np.ceil(7500000 // self.skipping // self.batch_size // self.slice_length)\n\n        stride = 1\n\n        # We use tensorflow datasets as it is a very efficient way to load data when using big data chunks\n        # Step 1: Create a dataset from the CSV file\n        # Read data, create windows of equal length and batch them into chunks\n        dataset1 = tf.data.experimental.CsvDataset(\n            filenames=self.path_area1,\n            record_defaults=column_defaults,\n            header=True,\n            exclude_cols=[0]\n        )\n        dataset1 = dataset1.map(lambda *cols: tf.stack(cols, axis=0))\n        dataset1 = dataset1.window(size=slice_length, shift=stride, drop_remainder=True)\n        dataset1 = dataset1.flat_map(lambda *window: tf.data.Dataset.zip(window).batch(slice_length))\n\n        dataset2 = tf.data.experimental.CsvDataset(\n            filenames=self.path_area2,\n            record_defaults=column_defaults,\n            header=True,\n            exclude_cols=[0]\n        )\n\n        dataset2 = dataset2.map(lambda *cols: tf.stack(cols, axis=0))\n        dataset2 = dataset2.window(size=slice_length, shift=stride, drop_remainder=True)\n        dataset2 = dataset2.flat_map(lambda *window: tf.data.Dataset.zip(window).batch(slice_length))\n\n        probabilities = [0.5, 0.5]\n\n        # Sample from the two datasets\n        sampled_dataset = tf.data.experimental.sample_from_datasets([dataset1, dataset2], probabilities)\n\n        # Assuming `dataset` is your TensorFlow dataset\n        sampled_dataset = sampled_dataset.map(preprocess_for_autoencoder)\n\n        # Repeat the dataset indefinitely\n        sampled_dataset = sampled_dataset.repeat()\n\n        # Batch and prefetch for performance\n        batched_dataset = sampled_dataset.batch(batch_size, drop_remainder=True).prefetch(tf.data.experimental.AUTOTUNE)\n        self.sampled_dataset = batched_dataset\n\n    def train(self,\n              latent_dim=24,\n              patience = 20,\n              max_epochs=100):\n        self.epochs = max_epochs\n        self.latent_dim = latent_dim\n        self.patience = patience\n\n        inputs = Input(shape=(self.slice_length,\n                              self.num_features))\n\n        # Encoder - dimensionality reduction\n        encoder_lstm1 = LSTM(128, return_sequences=True)(inputs)\n        encoder_lstm1 = Dropout(0.1)(encoder_lstm1)\n        encoder_lstm2 = LSTM(64, return_sequences=False)(encoder_lstm1)\n        dropout_enc1 = Dropout(0.1)(encoder_lstm2)\n\n        # Here we are introducing a variational part\n        z_mean = Dense(latent_dim)(dropout_enc1)\n        z_log_var = Dense(latent_dim)(dropout_enc1)\n\n        def sampling(args):\n            z_mean, z_log_var = args\n            epsilon = K.random_normal(shape=(K.shape(z_mean)[0], latent_dim),\n                                      mean=0., stddev=1.)\n            return z_mean + K.exp(z_log_var) * epsilon\n\n        z = Lambda(sampling, output_shape=(latent_dim,))([z_mean, z_log_var])\n\n        # Decoder - reconstruct the data\n        repeat_z = RepeatVector(self.slice_length)(z)\n        decoder_lstm = LSTM(64, return_sequences=True)(repeat_z)\n        decoder_conv1 = Dropout(0.05)(decoder_lstm)\n        decoder_lstm2 = LSTM(128, return_sequences=True)(decoder_conv1)\n        outputs = TimeDistributed(Dense(self.num_features))(decoder_lstm2)\n\n        def vae_loss(x, x_decoded_mean):\n            xent_loss = K.mean(tf.keras.losses.mse(x, x_decoded_mean), axis=-1)\n            kl_loss = - 0.5 * K.mean(1 + z_log_var - K.square(z_mean) - K.exp(z_log_var), axis=-1)\n            return xent_loss + kl_loss\n\n        # note that we take them as a list here\n        self.model = Model(inputs, outputs)\n\n        self.model.add_loss(vae_loss(inputs, outputs))\n\n        self.model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.001))\n\n        # model callback local\n        model_checkpoint_callback = tf.keras.callbacks.ModelCheckpoint(\n            filepath=f\"vae_model.\"+'{epoch:02d}-{loss:.2f}.hdf5',  # Path where to save the model\n            save_weights_only=False,  # Set to `True` to save only weights, `False` to save the whole model\n            monitor='loss',  # Metric to monitor\n            mode='min',  # The model is saved when the monitored metric has minimized\n            save_best_only=False,  # Only the best model is saved\n            verbose=2  # Verbosity mode\n        )\n\n        # Fit the defined Neural Network from above\n        self.model.fit(self.sampled_dataset,\n                epochs=self.epochs,\n                steps_per_epoch=self.steps_per_epoch,\n                callbacks=[tf.keras.callbacks.EarlyStopping(monitor=\"loss\", patience=self.patience, mode=\"min\"),\n                          model_checkpoint_callback\n                          ],\n                       verbose=2\n                )\n\n        # save the model\n        self.model.save(f\"vae_model.keras\")\n\n    def predict(self,\n                test_path,\n                submission_path):\n        self.test_path = test_path\n\n        df_test = pd.read_csv(test_path)\n\n        # split test data into controlAreas\n        segment_ids_1 = df_test[df_test['controlArea']==1]['test_data_segment_id'].unique()\n        segment_ids_2 = df_test[df_test['controlArea']==2]['test_data_segment_id'].unique()\n\n        test_predictions_1 =  pd.DataFrame( columns = ['anomaly'])\n        \n        # predict data per segment and calculate the anomaly score (MSE) for the reconstructed data\n        for segment_id in segment_ids_1:\n            print(\"predict segment \", segment_id)\n            start = time.time()\n            df_test_sample = df_test[df_test['test_data_segment_id']== segment_id].drop(columns=[ 'Datum_Uhrzeit_CET', 'test_data_segment_id'])\n\n            df_reconstructed_mean, _ = reconstruct_data_segment(df = df_test_sample, window_size = self.slice_length, model = self.model)\n\n            segment_pred = recontruction_to_anomalies(df_test_sample, df_reconstructed_mean)\n\n            test_predictions_1 = pd.concat([test_predictions_1, segment_pred])\n            end = time.time()\n            print(f\"finished segment {segment_id} in {end - start} s\")\n\n        test_predictions_1.to_csv(\"test_pred_1.csv\")\n\n        test_predictions_2 =  pd.DataFrame( columns = ['anomaly'])\n\n        for segment_id in segment_ids_2:\n            print(\"predict segment \", segment_id)\n            start = time.time()\n            df_test_sample = df_test[df_test['test_data_segment_id']== segment_id].drop(columns=[ 'Datum_Uhrzeit_CET', 'test_data_segment_id'])\n\n            df_reconstructed_mean, _ = reconstruct_data_segment(df = df_test_sample, window_size = self.slice_length, model = self.model)\n\n            segment_pred = recontruction_to_anomalies(df_test_sample, df_reconstructed_mean)\n\n            test_predictions_2 = pd.concat([test_predictions_2, segment_pred])\n            end = time.time()\n            print(f\"finished segment {segment_id} in {end - start} s\")\n        test_predictions_2.to_csv(\"test_pred_2.csv\")\n\n\n        submission = pd.concat([test_predictions_1, test_predictions_2])\n        submission.index = submission.index.set_names([\"id\"])\n        submission = submission.reset_index()\n\n        submission.to_csv(\"submission_vae.csv\", index=False)\n        submission.to_csv(submission_path, index=False)\n\n        return submission\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Traning the model on the train set","metadata":{}},{"cell_type":"code","source":"model = TimeSeries_LSTM_Variational_Autoencoder()\nmodel.configure_for_training(\n    path_area1 = \"final_area1_train.csv\",\n    path_area2 = \"final_area2_train.csv\",\n    skipping = 4,\n    slice_length = 32,\n    batch_size = 64,\n)\n\n# Note: The train time is quite long, thus we trained our models in a distributed worker strategy\n# with google colab TPUs, which requires to store pre-processed traning data on a google cloud storage bucket\nmodel.train(latent_dim=16,\n            patience = 10,\n            max_epochs = 100)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Applying the model on the test set","metadata":{}},{"cell_type":"code","source":"# Create submission based on ML model only\nsubmission = model.predict(test_path = \"final_test.csv\",\n                           submission_path = \"raw_submission.csv\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Afterprocessing the submission data","metadata":{}},{"cell_type":"markdown","source":"We tried to clean the output data to make the results and plots more reasonable, as we noticed that the models have some issues that are technically not correct but can be minimized by rule based afterprocessing.","metadata":{}},{"cell_type":"code","source":"# Filter out anomaly windows that are too short and would not make sense in a productive environment.\n\ndef filter_time_series(df2, column, min_range_length):\n    df = df2.copy()\n    ranges = []\n    current_range = []\n    # Iterate through column 'A'\n    for index in df.index:\n        if df.loc[index,column] == 1:\n            current_range.append(index)\n        else:\n            if len(current_range) > 0:\n                ranges.append(current_range)\n            current_range = []\n \n    filtered_series = []\n    for range in ranges:\n        if len(range) >= min_range_length:\n            filtered_series.extend(range)\n        # else:\n        #     filtered_series.extend([0] * len(range))\n\n    df.loc[:,\"old_\" + column] = df.loc[:, column].copy()\n    df.loc[:, column] = 0\n    df.loc[filtered_series, column] = 1\n    return df\n\n# Close gaps between big blocks of anomalies\ndef close_gaps(df2, max_thr):\n    df = df2.copy()\n    ranges = []\n    current_range = []\n    # Iterate through column 'A'\n    for index in df.index:\n        if df.loc[index,\"anomaly\"] == 1:\n            current_range.append(index)\n        else:\n            if len(current_range) > 0:\n                ranges.append(current_range)\n            current_range = []\n\n\n    new_windows = []\n    for index, current_range in enumerate(ranges):\n        if index == len(ranges)-1:\n            continue\n        else:\n            len1 = len(current_range)\n            len2 = len(ranges[index+1])\n            max_ind_1 = current_range[-1]\n            min_ind_2 = ranges[index+1][0]\n            len_gap = min_ind_2 - max_ind_1\n            if (len_gap < len1) & (len_gap < len2) & (len_gap < max_thr):\n                closed_window = range(max_ind_1 + 1, min_ind_2)\n                new_windows.extend(closed_window)\n\n    # print(f\"added: {(df.loc[new_windows,'anomaly'] == 0).sum()}\")\n    df.loc[new_windows, \"anomaly\"] = 1\n    return df\n\ndef afterprocessing(df, threshold):\n\n    submission = df.copy()\n    # Smooth window - count the anomalies around. If there are enough \n    # and the anomaly score is big enough compared to the threshold, mark anomaly\n    submission[\"count_10_around\"] = submission[\"anomaly\"].rolling(window = 10, min_periods = 1, center=True).sum().astype(int)\n    submission[\"anomaly_1\"] = (submission[\"anomaly\"] == 1) | ( (submission[\"score\"] > 0.9 *threshold) & (submission[\"count_10_around\"]>=2) ).astype(int) | (submission[\"count_10_around\"]>=8)\n    submission[\"count_10_around_1\"] = submission[\"anomaly_1\"].rolling(window = 10, min_periods = 1, center=True).sum().astype(int)\n    submission[\"anomaly_2\"] = (submission[\"anomaly_1\"] == 1) | ( (submission[\"score\"] > 0.95 * threshold) & (submission[\"count_10_around_1\"]>=7) ).astype(int) | (submission[\"count_10_around\"]>=9)\n\n    # Delete singles - Anomaly windows with less than 4 indices don't make sense\n    submission[\"anomaly_single\"] = (submission[\"count_10_around\"] < 4) & (submission[\"anomaly\"] == 1).astype(int)\n\n    submission[\"anomaly_old\"] = submission[\"anomaly\"]\n\n    # combine the found anomalies\n    submission[\"anomaly\"] = (submission[\"anomaly\"] | submission[\"anomaly_1\"] | submission[\"anomaly_2\"]).astype(int) - submission[\"anomaly_single\"]\n\n    return submission[[\"id\",\"anomaly\",\"score\"]]\n\ndef clean_submission(submission, threshold):\n    #print(f\"before: {(submission['anomaly'] > threshold).sum()}\")\n    if \"score\" not in submission.columns:\n      submission[\"score\"] = submission[\"anomaly\"].copy()\n      submission[\"anomaly\"] = (submission[\"score\"]>threshold).astype(int)\n    print(f\"before: {submission['anomaly'].sum()}\")\n    for segment in X_test[\"test_data_segment_id\"].unique():\n        submission.loc[X_test[\"test_data_segment_id\"]==segment,:] = afterprocessing(submission.loc[X_test[\"test_data_segment_id\"]==segment,:], threshold)\n        submission.loc[X_test[\"test_data_segment_id\"]==segment,:] = afterprocessing(submission.loc[X_test[\"test_data_segment_id\"]==segment,:], threshold)\n        submission.loc[X_test[\"test_data_segment_id\"]==segment,:] = filter_time_series(submission.loc[X_test[\"test_data_segment_id\"]==segment,:],\"anomaly\",20)\n        submission.loc[X_test[\"test_data_segment_id\"]==segment,:] = close_gaps(submission.loc[X_test[\"test_data_segment_id\"]==segment,:],300)\n    submission.loc[(X_test_filter[\"participationCMO\"]==0) | (X_test_filter[\"participationIN\"]==0),\"anomaly\"] = 0\n    print(f\"after: {submission['anomaly'].sum()}\")\n    return(submission)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission1 = pd.read_csv(\"raw_submission.csv\") \nthreshold = submission1[\"anomaly\"].quantile(0.75) # choose a reasonable quantile to set the threshold - found by testing and plotting\nsubmission2 = clean_submission(submission1, threshold)\nsubmission3 = clean_submission(submission2, threshold).set_index('id')\nsubmission = submission3.copy()","metadata":{"execution":{"iopub.status.busy":"2024-02-13T17:34:10.136764Z","iopub.execute_input":"2024-02-13T17:34:10.137620Z","iopub.status.idle":"2024-02-13T17:34:10.472168Z","shell.execute_reply.started":"2024-02-13T17:34:10.137585Z","shell.execute_reply":"2024-02-13T17:34:10.470908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Combine ML results with from the rule based approach\nWe found that only the rule for \"(near) constant values\" needs to be added to the output of cleaned ML-based prediciton, the other rule-based anomalies are already captured quite-well by integration them as featues in the ML model.","metadata":{}},{"cell_type":"code","source":"submission[\"anomaly\"] = (submission[\"anomaly\"].sort_index() | rb_test[[\"id\", \"flagNearConstantValue\"]].set_index('id').sort_index()[\"flagNearConstantValue\"]).astype(int) #(submission[\"anomaly\"] | rb_test[\"flagNearConstantValue\"]).astype(int)\nsubmission[\"anomaly\"].to_csv(\"submission.csv\", index=True)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot the results for visual inspection of test prediciton","metadata":{}},{"cell_type":"code","source":"# Put in the submission file and plot the anomaly windows together with selected columns of the test dataframe\ndef plot_anomaly_window(df,cols,anomaly):\n    plt.figure(figsize=(25,10))\n\n    # Plot the time series\n    for col in cols:\n        plt.plot(df.index, df[col], label=col)\n\n    window_starts = []\n    window_ends = []\n    col_count=0\n\n    in_window = False\n\n    for index in anomaly.index:\n        value = anomaly.loc[index,\"anomaly\"]\n        if value == 1 and not in_window:\n            window_starts.append(index)\n            in_window = True\n        elif value == 0 and in_window:\n            window_ends.append(index - 1)\n            in_window = False\n\n    # If the last window continues until the end of the column, add it\n    if in_window:\n        window_ends.append(anomaly.index[-1] - 1)\n\n    if(len(window_starts)>0):\n      for start, end in zip(window_starts, window_ends):\n            plt.axvspan(start, end, facecolor='g', alpha=0.3)\n\n    # Set labels and title\n    plt.xlabel('Date')\n    plt.ylabel('Value')\n    plt.title('Time Series with Flagged Intervals')\n\n    # Add legend\n    plt.legend()\n\n    # Show the plot\n    plt.show()\n\ndef plot_test_anomalies(X_test, submission, columns=[\"aFRRactivation\",\"aFRRrequest\",\"correctedDemand\"]):\n    for segment in X_test[\"test_data_segment_id\"].unique():\n        print(segment)\n        segment_index = X_test[X_test[\"test_data_segment_id\"]==segment].index\n        plot_anomaly_window(X_test.iloc[segment_index,:],columns,submission.iloc[segment_index,:])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Columns of interest can be adjusted here\nplot_test_anomalies(X_test, submission, columns=[\"aFRRactivation\",\"aFRRrequest\",\"correctedDemand\"])","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}