{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"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":"none","dataSources":[{"sourceId":56537,"databundleVersionId":8015876,"sourceType":"competition"},{"sourceId":177694049,"sourceType":"kernelVersion"}],"dockerImageVersionId":30698,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Autoencoder as a Dimension Reduction technique\n\nAlong with PCA the embedding of the initial features (vector representation in low dimensional space) could be used.\n","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport tensorflow as tf\nfrom pathlib import Path","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-05-18T15:07:39.571885Z","iopub.execute_input":"2024-05-18T15:07:39.572328Z","iopub.status.idle":"2024-05-18T15:07:56.995403Z","shell.execute_reply.started":"2024-05-18T15:07:39.572294Z","shell.execute_reply":"2024-05-18T15:07:56.993679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fix the seed","metadata":{}},{"cell_type":"code","source":"import os \nimport random\nimport numpy as np \n\nDEFAULT_RANDOM_SEED = 1000\n\nrandom.seed(DEFAULT_RANDOM_SEED)\nos.environ['PYTHONHASHSEED'] = str(DEFAULT_RANDOM_SEED)\nnp.random.seed(DEFAULT_RANDOM_SEED)\ntf.random.set_seed(DEFAULT_RANDOM_SEED)","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:07:56.997742Z","iopub.execute_input":"2024-05-18T15:07:56.998454Z","iopub.status.idle":"2024-05-18T15:07:57.007273Z","shell.execute_reply.started":"2024-05-18T15:07:56.998418Z","shell.execute_reply":"2024-05-18T15:07:57.006006Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading data\n\nTo train Variational Autoencoder the random subsample of the all data will be used. It's about 1M rows.\n","metadata":{}},{"cell_type":"code","source":"data_dir = Path(\"/kaggle/input/create-random-sample\")\n\ndf = pd.concat(\n    pd.read_parquet(parquet_file)\n    for i,parquet_file in enumerate(data_dir.glob('*.parquet'))\n)","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:07:57.055333Z","iopub.execute_input":"2024-05-18T15:07:57.05615Z","iopub.status.idle":"2024-05-18T15:08:47.662016Z","shell.execute_reply.started":"2024-05-18T15:07:57.056105Z","shell.execute_reply":"2024-05-18T15:08:47.660874Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features = [\n    \"state_t\",\n    \"state_q0001\",\n    \"state_q0002\",\n    \"state_q0003\",\n    \"state_u\",\n    \"state_v\",\n    \"state_ps\",\n    \"pbuf_SOLIN\",\n    \"pbuf_LHFLX\",\n    \"pbuf_SHFLX\",\n    \"pbuf_TAUX\",\n    \"pbuf_TAUY\",\n    \"pbuf_COSZRS\",\n    \"cam_in_ALDIF\",\n    \"cam_in_ALDIR\",\n    \"cam_in_ASDIF\",\n    \"cam_in_ASDIR\",\n    \"cam_in_LWUP\",\n    \"cam_in_ICEFRAC\",\n    \"cam_in_LANDFRAC\",\n    \"cam_in_OCNFRAC\",\n    \"cam_in_SNOWHLAND\",\n    \"pbuf_ozone\",\n#     \"pbuf_CH4\",\n#     \"pbuf_N2O\",\n]\n\ntargets = [\n    \"ptend_t\",\n    \"ptend_q0001\",\n    \"ptend_q0002\",\n    \"ptend_q0003\",\n    \"ptend_u\",\n    \"ptend_v\",\n    \"cam_out_NETSW\",\n    \"cam_out_FLWDS\",\n    \"cam_out_PRECSC\",\n    \"cam_out_PRECC\",\n    \"cam_out_SOLS\",\n    \"cam_out_SOLL\",\n    \"cam_out_SOLSD\",\n    \"cam_out_SOLLD\",\n]","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:08:47.663301Z","iopub.execute_input":"2024-05-18T15:08:47.663627Z","iopub.status.idle":"2024-05-18T15:08:47.671574Z","shell.execute_reply.started":"2024-05-18T15:08:47.663599Z","shell.execute_reply":"2024-05-18T15:08:47.670313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# collect all feature and target column names with indices\nall_targets = []\nfor t in targets:\n    all_targets += list(filter(lambda x: x.startswith(t), df.columns))\n\nall_targets = list(set(all_targets))\n\nall_features = []\nfor f in features:\n    all_features += list(filter(lambda x: x.startswith(f), df.columns))","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:08:47.673078Z","iopub.execute_input":"2024-05-18T15:08:47.673513Z","iopub.status.idle":"2024-05-18T15:08:47.698211Z","shell.execute_reply.started":"2024-05-18T15:08:47.673484Z","shell.execute_reply":"2024-05-18T15:08:47.696911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Autoencoder model\n\nRead more about autoencoders [here](https://www.tensorflow.org/tutorials/generative/autoencoder).\n\nUsually the embedding approach is used to represent complex data such as pictures, text etc. In this notebook we will build shorter representation of longer vector.\n","metadata":{}},{"cell_type":"code","source":"class Autoencoder(tf.keras.models.Model):\n    def __init__(self, input_dim: int, latent_dim: int):\n        \n        # assure that encoder and decoder layers size change monotonically\n        assert input_dim // 2 > latent_dim\n        \n        super(Autoencoder, self).__init__()\n\n        self.encoder = tf.keras.Sequential([\n            tf.keras.layers.Dense(input_dim, activation=\"silu\"),\n            tf.keras.layers.Dense(input_dim // 2, activation=\"silu\"),\n            tf.keras.layers.Dense(input_dim // 3, activation=\"silu\"),\n            tf.keras.layers.Dense(latent_dim, activation=\"silu\")])\n\n        self.decoder = tf.keras.Sequential([\n            tf.keras.layers.Dense(input_dim // 2, activation=\"silu\"),\n            tf.keras.layers.Dense(input_dim // 3, activation=\"silu\"),\n            tf.keras.layers.Dense(input_dim, activation=\"relu\"),\n        ])\n\n    def call(self, x):\n        encoded = self.encoder(x)\n        decoded = self.decoder(encoded)\n        return decoded","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:08:47.699559Z","iopub.execute_input":"2024-05-18T15:08:47.700052Z","iopub.status.idle":"2024-05-18T15:08:47.712192Z","shell.execute_reply.started":"2024-05-18T15:08:47.700004Z","shell.execute_reply":"2024-05-18T15:08:47.711017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Preparation\n\nBefore being used as an input to neural net the data should be scaled and aometimes normalized. Let's apply simplified min-max scaling and filter out outliers using z-score.","metadata":{}},{"cell_type":"code","source":"def scale_minmax(df: pd.DataFrame, mins: pd.Series = None, maxs: pd.Series = None) -> tuple:\n    mins = mins if mins else df.min() \n    maxs = maxs if maxs else df.max()\n\n    scaled = (df - mins) / (maxs - mins)\n    \n    return scaled, mins, maxs\n\ndef scale_minmax_inverse(df: pd.DataFrame, mins: pd.Series, maxs: pd.Series) -> pd.DataFrame:\n    return df * (maxs - mins) + mins\n\ndef normalize_zscore(df: pd.DataFrame, means: pd.Series=None, stds:pd.Series=None) -> tuple:\n    means = means if means else df.mean()\n    stds = stds if stds else df.std()\n    \n    normalized = (df - means) / stds\n    \n    return normalized, means, stds\n\ndef normalize_zscore_inverse(df: pd.DataFrame, means:pd.Series, maxs: pd.Series) -> pd.DataFrame:\n    return df * std + means","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:08:47.713438Z","iopub.execute_input":"2024-05-18T15:08:47.713789Z","iopub.status.idle":"2024-05-18T15:08:47.734998Z","shell.execute_reply.started":"2024-05-18T15:08:47.713761Z","shell.execute_reply":"2024-05-18T15:08:47.733857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaled, mins, maxs = scale_minmax(df[all_features])\nnormalized, means, stds = normalize_zscore(scaled)\n\nfiltered = scaled[(normalized.abs() < 4).all(axis=1)]","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:08:47.736743Z","iopub.execute_input":"2024-05-18T15:08:47.737222Z","iopub.status.idle":"2024-05-18T15:09:02.632403Z","shell.execute_reply.started":"2024-05-18T15:08:47.737181Z","shell.execute_reply":"2024-05-18T15:09:02.631227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# shuffle data\nfiltered = filtered.sample(frac=1, random_state=42)","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:09:02.637747Z","iopub.execute_input":"2024-05-18T15:09:02.638157Z","iopub.status.idle":"2024-05-18T15:09:05.64337Z","shell.execute_reply.started":"2024-05-18T15:09:02.638125Z","shell.execute_reply":"2024-05-18T15:09:05.642239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# observe data that actually comes to our model\nfiltered","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:09:05.644889Z","iopub.execute_input":"2024-05-18T15:09:05.645245Z","iopub.status.idle":"2024-05-18T15:09:05.771598Z","shell.execute_reply.started":"2024-05-18T15:09:05.645216Z","shell.execute_reply":"2024-05-18T15:09:05.77038Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Features\n\nAs an example we will use only subset of all available features.","metadata":{}},{"cell_type":"code","source":"states_t = list(filter(lambda x: x.startswith(\"state_t\"), all_features))","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:09:05.772796Z","iopub.execute_input":"2024-05-18T15:09:05.773149Z","iopub.status.idle":"2024-05-18T15:09:05.779404Z","shell.execute_reply.started":"2024-05-18T15:09:05.773121Z","shell.execute_reply":"2024-05-18T15:09:05.778111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train-val split\n\nSplit data using 80/20 rule","metadata":{}},{"cell_type":"code","source":"TRAIN_SIZE = int(filtered.shape[0] * 0.8)\nBATCH_SIZE = 1024","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:09:05.780989Z","iopub.execute_input":"2024-05-18T15:09:05.781463Z","iopub.status.idle":"2024-05-18T15:09:05.794349Z","shell.execute_reply.started":"2024-05-18T15:09:05.781421Z","shell.execute_reply":"2024-05-18T15:09:05.793111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Building tf.data Dataset","metadata":{}},{"cell_type":"code","source":"# in the simple version of Autoencoder it is required to duplicate inputs as outputs\n\ntrain_sample = filtered[states_t][:TRAIN_SIZE]\ntest_sample = filtered[states_t][TRAIN_SIZE:]\n\ntrain_dataset = tf.data.Dataset.from_tensor_slices((train_sample, train_sample)).batch(BATCH_SIZE).prefetch(5).cache()\ntest_dataset = tf.data.Dataset.from_tensor_slices((test_sample, test_sample)).batch(BATCH_SIZE).prefetch(5).cache()","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:09:05.795908Z","iopub.execute_input":"2024-05-18T15:09:05.796352Z","iopub.status.idle":"2024-05-18T15:09:07.022234Z","shell.execute_reply.started":"2024-05-18T15:09:05.796321Z","shell.execute_reply":"2024-05-18T15:09:07.021064Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Initializing the model\n\nLet's initialize model with default parameters.\n\nThe model will optimize mean absolute difference between input sequence and reconstructed one. We will also track percentage error MAPE as a main metric.","metadata":{}},{"cell_type":"code","source":"autoencoder = Autoencoder(input_dim=60, latent_dim=16)\n\noptimizer = tf.keras.optimizers.Adam(1e-4)\nloss = tf.keras.losses.MeanAbsoluteError()\nautoencoder.compile(optimizer=optimizer, loss=loss, metrics=['mape'])","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:09:07.024042Z","iopub.execute_input":"2024-05-18T15:09:07.024409Z","iopub.status.idle":"2024-05-18T15:09:07.061245Z","shell.execute_reply.started":"2024-05-18T15:09:07.024379Z","shell.execute_reply":"2024-05-18T15:09:07.05993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training the model\n\nThe current light-weight Autoencoder implementation doesn't requre GPU to train.","metadata":{}},{"cell_type":"code","source":"history = autoencoder.fit(\n    x=train_dataset,\n    epochs=200,\n    validation_data=test_dataset,\n    shuffle=True,\n    callbacks=[\n        tf.keras.callbacks.EarlyStopping(\n            monitor=\"val_mape\",\n            min_delta=0.0,\n            patience=10\n        ),\n        tf.keras.callbacks.ModelCheckpoint(\n            \"/kaggle/working/checkpoint.model.keras\",\n            monitor=\"val_total_loss\",\n            verbose=0,\n            save_best_only=True,\n            save_weights_only=False,\n            mode=\"min\",\n            save_freq=\"epoch\",\n            initial_value_threshold=None,\n        )\n    ],\n    verbose=1\n)","metadata":{"execution":{"iopub.status.busy":"2024-05-18T15:09:07.062306Z","iopub.execute_input":"2024-05-18T15:09:07.062636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Collect predictions on validation set","metadata":{}},{"cell_type":"code","source":"embeddings_val = autoencoder.encoder(test_sample)\npreds = autoencoder.decoder(embeddings_val)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizing the results","metadata":{}},{"cell_type":"code","source":"import plotly.express as px\nimport plotly.graph_objects as go","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Learning curve","metadata":{}},{"cell_type":"code","source":"history_to_vis = history.history.copy()\n\nhistory_to_vis['mape'] = np.array(history_to_vis['mape']) / 100\nhistory_to_vis['val_mape'] = np.array(history_to_vis['val_mape']) / 100","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.line(history_to_vis, title='Learning curve')\nfig.update_xaxes(title={\"text\": \"Epoch\"})","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Embeddings","metadata":{}},{"cell_type":"code","source":"px.line(embeddings_val[1000:1020].numpy().T)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizing Ground truths against predictions","metadata":{}},{"cell_type":"markdown","source":"## Comparison graph","metadata":{}},{"cell_type":"code","source":"fig = go.Figure()\n\nbatch_number = 50\nitem_number = 100\n\nfig.add_trace(go.Scatter(y=preds[BATCH_SIZE * batch_number + item_number].numpy().T, mode='lines', name=\"pred\"))\nfig.add_trace(go.Scatter(y=test_sample.iloc[BATCH_SIZE * batch_number + item_number], mode='lines', name='true'))\nfig.update_layout(title=dict(text=\"Reconstructed input vs ground truth\"))\n\nfig.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Error analysis","metadata":{}},{"cell_type":"code","source":"fig = px.box((preds - test_sample)[:10000])\nfig.update_layout(title=dict(text=\"Errors distribution per feature\"))","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}