{
 "cells": [
  {
   "cell_type": "raw",
   "id": "5f5c077c",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/adriangb/scikeras/blob/docs-deploy/refs/heads/master/notebooks/AutoEncoders.ipynb\"><img src=\"https://www.tensorflow.org/images/colab_logo_32px.png\">Run in Google Colab</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76a4d363",
   "metadata": {},
   "source": [
    "# Autoencoders in SciKeras\n",
    "\n",
    "Autencoders are an approach to use nearual networks to distill data into it's most important features, thereby compressing the data.\n",
    "We will be following the [Keras tutorial](https://blog.keras.io/building-autoencoders-in-keras.html) on the topic, which goes much more in depth and breadth than we will here.\n",
    "You are highly encouraged to check out that tutorial if you want to learn about autoencoders in the general sense.\n",
    "\n",
    "## Table of contents\n",
    "\n",
    "* [1. Setup](#1.-Setup)\n",
    "* [2. Data](#2.-Data)\n",
    "* [3. Define Keras Model](#3.-Define-Keras-Model)\n",
    "* [4. Training](#4.-Training)\n",
    "* [5. Explore Results](#5.-Explore-Results)\n",
    "* [6. Deep AutoEncoder](#6.-Deep-AutoEncoder)\n",
    "\n",
    "## 1. Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a16b3c5c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:48:36.379232Z",
     "iopub.status.busy": "2021-11-29T20:48:36.372314Z",
     "iopub.status.idle": "2021-11-29T20:48:39.167543Z",
     "shell.execute_reply": "2021-11-29T20:48:39.166367Z"
    }
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import scikeras\n",
    "except ImportError:\n",
    "    !python -m pip install scikeras"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "410d8f2a",
   "metadata": {},
   "source": [
    "Silence TensorFlow logging to keep output succinct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "47c66c84",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:48:39.174034Z",
     "iopub.status.busy": "2021-11-29T20:48:39.170690Z",
     "iopub.status.idle": "2021-11-29T20:48:39.177530Z",
     "shell.execute_reply": "2021-11-29T20:48:39.178361Z"
    }
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "from tensorflow import get_logger\n",
    "get_logger().setLevel('ERROR')\n",
    "warnings.filterwarnings(\"ignore\", message=\"Setting the random state for TF\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "880a2f08",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:48:39.188429Z",
     "iopub.status.busy": "2021-11-29T20:48:39.182026Z",
     "iopub.status.idle": "2021-11-29T20:48:39.656703Z",
     "shell.execute_reply": "2021-11-29T20:48:39.657252Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scikeras.wrappers import KerasClassifier, KerasRegressor\n",
    "from tensorflow import keras"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34724672",
   "metadata": {},
   "source": [
    "## 2. Data\n",
    "\n",
    "We load the dataset from the Keras tutorial. The dataset consists of images of cats and dogs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "215f0b10",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:48:39.665210Z",
     "iopub.status.busy": "2021-11-29T20:48:39.663308Z",
     "iopub.status.idle": "2021-11-29T20:48:40.328300Z",
     "shell.execute_reply": "2021-11-29T20:48:40.329141Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(60000, 784)\n",
      "(10000, 784)\n"
     ]
    }
   ],
   "source": [
    "from tensorflow.keras.datasets import mnist\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "(x_train, _), (x_test, _) = mnist.load_data()\n",
    "x_train = x_train.astype('float32') / 255.\n",
    "x_test = x_test.astype('float32') / 255.\n",
    "x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:])))\n",
    "x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:])))\n",
    "print(x_train.shape)\n",
    "print(x_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "718a70f1",
   "metadata": {},
   "source": [
    "## 3. Define Keras Model\n",
    "\n",
    "We will be defining a very simple autencoder. We define _three_ model architectures:\n",
    "\n",
    "1. An encoder: a series of densly connected layers culminating in an \"output\" layer that determines the encoding dimensions.\n",
    "2. A decoder: takes the output of the encoder as it's input and reconstructs the original data.\n",
    "3. An autoencoder: a chain of the encoder and decoder that directly connects them for training purposes.\n",
    "\n",
    "The only variable we give our model is the encoding dimensions, which will be a hyperparemter of our final transformer.\n",
    "\n",
    "The encoder and decoder are views to the first/last layers of the autoencoder model.\n",
    "They'll be directly used in `transform` and `inverse_transform`, so we'll create some SciKeras models with those layers\n",
    "and save them as in `encoder_model_` and `decoder_model_`. All three models are created within `_keras_build_fn`.\n",
    "\n",
    "For a background on chaining Functional Models like this, see [All models are callable](https://keras.io/guides/functional_api/#all-models-are-callable-just-like-layers) in the Keras docs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0ef8e36a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:48:40.344847Z",
     "iopub.status.busy": "2021-11-29T20:48:40.332977Z",
     "iopub.status.idle": "2021-11-29T20:48:40.361381Z",
     "shell.execute_reply": "2021-11-29T20:48:40.360300Z"
    }
   },
   "outputs": [],
   "source": [
    "from typing import Dict, Any\n",
    "\n",
    "from sklearn.base import TransformerMixin\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from scikeras.wrappers import BaseWrapper\n",
    "\n",
    "\n",
    "class AutoEncoder(BaseWrapper, TransformerMixin):\n",
    "    \"\"\"A class that enables transform and fit_transform.\n",
    "    \"\"\"\n",
    "\n",
    "    encoder_model_: BaseWrapper\n",
    "    decoder_model_: BaseWrapper\n",
    "    \n",
    "    def _keras_build_fn(self, encoding_dim: int, meta: Dict[str, Any]):\n",
    "        n_features_in = meta[\"n_features_in_\"]\n",
    "\n",
    "        encoder_input = keras.Input(shape=(n_features_in,))\n",
    "        encoder_output = keras.layers.Dense(encoding_dim, activation='relu')(encoder_input)\n",
    "        encoder_model = keras.Model(encoder_input, encoder_output)\n",
    "\n",
    "        decoder_input = keras.Input(shape=(encoding_dim,))\n",
    "        decoder_output = keras.layers.Dense(n_features_in, activation='sigmoid', name=\"decoder\")(decoder_input)\n",
    "        decoder_model = keras.Model(decoder_input, decoder_output)\n",
    "        \n",
    "        autoencoder_input = keras.Input(shape=(n_features_in,))\n",
    "        encoded_img = encoder_model(autoencoder_input)\n",
    "        reconstructed_img = decoder_model(encoded_img)\n",
    "\n",
    "        autoencoder_model = keras.Model(autoencoder_input, reconstructed_img)\n",
    "\n",
    "        self.encoder_model_ = BaseWrapper(encoder_model, verbose=self.verbose)\n",
    "        self.decoder_model_ = BaseWrapper(decoder_model, verbose=self.verbose)\n",
    "\n",
    "        return autoencoder_model\n",
    "    \n",
    "    def _initialize(self, X, y=None):\n",
    "        X, _ = super()._initialize(X=X, y=y)\n",
    "        # since encoder_model_ and decoder_model_ share layers (and their weights)\n",
    "        # X_tf here come from random weights, but we only use it to initialize our models\n",
    "        X_tf = self.encoder_model_.initialize(X).predict(X)\n",
    "        self.decoder_model_.initialize(X_tf)\n",
    "        return X, X\n",
    "\n",
    "    def initialize(self, X):\n",
    "        self._initialize(X=X, y=X)\n",
    "        return self\n",
    "\n",
    "    def fit(self, X, *, sample_weight=None) -> \"AutoEncoder\":\n",
    "        super().fit(X=X, y=X, sample_weight=sample_weight)\n",
    "        # at this point, encoder_model_ and decoder_model_\n",
    "        # are both \"fitted\" because they share layers w/ model_\n",
    "        # which is fit in the above call\n",
    "        return self\n",
    "\n",
    "    def score(self, X) -> float:\n",
    "        # Note: we use 1-MSE as the score\n",
    "        # With MSE, \"larger is better\", but Scikit-Learn\n",
    "        # always maximizes the score (e.g. in GridSearch)\n",
    "        return 1 - mean_squared_error(self.predict(X), X)\n",
    "\n",
    "    def transform(self, X) -> np.ndarray:\n",
    "        X: np.ndarray = self.feature_encoder_.transform(X)\n",
    "        return self.encoder_model_.predict(X)\n",
    "\n",
    "    def inverse_transform(self, X_tf: np.ndarray):\n",
    "        X: np.ndarray = self.decoder_model_.predict(X_tf)\n",
    "        return self.feature_encoder_.inverse_transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bd91961",
   "metadata": {},
   "source": [
    "Next, we wrap the Keras Model with Scikeras. Note that for our encoder/decoder estimators, we do not need to provide a loss function since no training will be done.\n",
    "We do however need to have the `fit_model` and `encoding_dim` so that these will be settable by `BaseWrapper.set_params`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d56911a1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:48:40.366818Z",
     "iopub.status.busy": "2021-11-29T20:48:40.366189Z",
     "iopub.status.idle": "2021-11-29T20:48:40.369784Z",
     "shell.execute_reply": "2021-11-29T20:48:40.370822Z"
    }
   },
   "outputs": [],
   "source": [
    "autoencoder = AutoEncoder(\n",
    "    loss=\"binary_crossentropy\",\n",
    "    encoding_dim=32,\n",
    "    random_state=0,\n",
    "    epochs=5,\n",
    "    verbose=False,\n",
    "    optimizer=\"adam\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ac5714e",
   "metadata": {},
   "source": [
    "## 4. Training\n",
    "\n",
    "To train the model, we pass the input images as both the features and the target.\n",
    "This will train the layers to compress the data as accurately as possible between the encoder and decoder.\n",
    "Note that we only pass the `X` parameter, since we defined the mapping `y=X` in `KerasTransformer.fit` above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8e63c86a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:48:40.375769Z",
     "iopub.status.busy": "2021-11-29T20:48:40.375133Z",
     "iopub.status.idle": "2021-11-29T20:49:29.061017Z",
     "shell.execute_reply": "2021-11-29T20:49:29.061547Z"
    }
   },
   "outputs": [],
   "source": [
    "_ = autoencoder.fit(X=x_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d6ed7f2",
   "metadata": {},
   "source": [
    "Next, we round trip the test dataset and explore the performance of the autoencoder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7ba6a6d0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:49:29.065242Z",
     "iopub.status.busy": "2021-11-29T20:49:29.064574Z",
     "iopub.status.idle": "2021-11-29T20:49:30.524722Z",
     "shell.execute_reply": "2021-11-29T20:49:30.525268Z"
    }
   },
   "outputs": [],
   "source": [
    "roundtrip_imgs = autoencoder.inverse_transform(autoencoder.transform(x_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c16d5fc",
   "metadata": {},
   "source": [
    "## 5. Explore Results\n",
    "\n",
    "Let's compare our inputs to lossy decoded outputs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bfd16bdb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:49:30.534820Z",
     "iopub.status.busy": "2021-11-29T20:49:30.533974Z",
     "iopub.status.idle": "2021-11-29T20:49:40.659342Z",
     "shell.execute_reply": "2021-11-29T20:49:40.658355Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Matplotlib is building the font cache; this may take a moment.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABG0AAADnCAYAAACkCqtqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABCP0lEQVR4nO3dedxVU///8RVCaVAphGYZGillikq3IcoUusuYeQy3ebiR6XuHCCkhQypTZUwUGUK3X0lpVqmkedAkRNfvj/vh471W19md63TOufZ1zuv512db6zpnt/dZ++yzrc/6lCooKHAAAAAAAACIl+2KewcAAAAAAACwJR7aAAAAAAAAxBAPbQAAAAAAAGKIhzYAAAAAAAAxxEMbAAAAAACAGOKhDQAAAAAAQAztUJTOpUqVoj54MSkoKCiVjtfhHBarFQUFBVXT8UKcx+LDWMwJjMUcwFjMCYzFHMBYzAmMxRzAWMwJhY5FZtoA2TO/uHcAgHOOsQjEBWMRiAfGIhAPhY5FHtoAAAAAAADEEA9tAAAAAAAAYoiHNgAAAAAAADHEQxsAAAAAAIAY4qENAAAAAABADPHQBgAAAAAAIIZ4aAMAAAAAABBDPLQBAAAAAACIoR2KeweQn2644QaLy5Qp47U1btzY4k6dOiV8jb59+1r81VdfeW0DBw7c1l0EAAAAAKBYMdMGAAAAAAAghnhoAwAAAAAAEEM8tAEAAAAAAIgh1rRB1rz66qsWR61VozZv3pyw7dJLL7W4Xbt2Xtunn35q8YIFC5LdRRSz+vXre9szZsywuHv37hY/8cQTWdunfLbLLrtY/NBDD1msY8855yZMmGDxGWec4bXNnz8/Q3sHAABQPCpVqmRxjRo1kvqb8J7ouuuus3jKlCkWz5o1y+s3adKkVHYROYSZNgAAAAAAADHEQxsAAAAAAIAYIj0KGaPpUM4lnxKlKTEffPCBxXXq1PH6dejQweK6det6bV27drX4wQcfTOp9UfwOOuggb1vT4xYuXJjt3cl7e+65p8UXX3yxxWHaYrNmzSw+6aSTvLY+ffpkaO+gDj74YIuHDRvmtdWqVStj73vsscd629OnT7f4xx9/zNj7Yuv0O9I5595++22Lr7rqKov79evn9fvzzz8zu2M5qFq1aha/9tprFn/55Zdev/79+1s8b968jO/XXypWrOhtH3XUURaPHDnS4k2bNmVtn4CS4MQTT7S4Y8eOXlvr1q0trlevXlKvF6Y91axZ0+Kddtop4d9tv/32Sb0+chczbQAAAAAAAGKIhzYAAAAAAAAxRHoU0qp58+YWn3rqqQn7TZ061eJwuuGKFSssXr9+vcU77rij12/cuHEWN2nSxGurUqVKknuMOGnatKm3vWHDBouHDx+e5b3JP1WrVvW2X3zxxWLaExTVcccdZ3HUFOt0C1NwunXrZnHnzp2zth/4H/3ue+qppxL2e/LJJy0eMGCA17Zx48b071iO0aoxzvn3NJqKtHTpUq9fcaVEaYU/5/xrvaa3zp49O/M7VsJUqFDB29aU+4YNG1ocVjEl1SzedFmFK6+80mJNBXfOuTJlylhcqlSpbX7fsEoqkCxm2gAAAAAAAMQQD20AAAAAAABiiIc2AAAAAAAAMVSsa9qEJaA1j3DRokVe26+//mrxoEGDLF6yZInXj3zc4qUlgsPcT8351vUXFi9enNRr/+tf//K2DzzwwIR933vvvaReE8VPc8K1DK1zzg0cODDbu5N3rrnmGotPOeUUr61FixZFfj0tJeucc9tt9/f/G5g0aZLFn332WZFfG74ddvj7K7x9+/bFsg/hWhnXX3+9xbvssovXpmtUITN0/O29994J+w0ZMsRivb9CYrvttpvFr776qtdWuXJli3UtoauvvjrzO5bAHXfcYXHt2rW9tksvvdRi7pu31LVrV4vvv/9+r22fffYp9G/CtW9WrlyZ/h1D2uj1sXv37hl9rxkzZlisv4WQPlpyXa/VzvlrrGqZduec27x5s8X9+vWz+IsvvvD6xeE6yUwbAAAAAACAGOKhDQAAAAAAQAwVa3pUz549ve1atWol9Xc6rXPdunVeWzannS1cuNDi8N8yfvz4rO1HnLzzzjsW61Q15/xztWrVqiK/dlg+tnTp0kV+DcTP/vvvb3GYThFOQUf6PfrooxbrNNFUnXbaaQm358+fb/FZZ53l9QvTbLB1bdq0sfiwww6zOPw+yqSw9LGmrZYtW9ZrIz0q/cLy7rfffntSf6eppwUFBWndp1x18MEHWxxOsVc9evTIwt5sqUGDBt62ppQPHz7ca+O7dUuaLvPYY49ZXKVKFa9fovHyxBNPeNua7p3KPS+SE6bCaKqTpriMHDnS6/fbb79ZvGbNGovD7ym9L/3www+9tilTplj83//+1+KJEyd6/TZu3Jjw9ZE8XU7BOX+M6b1m+JlIVsuWLS3+448/vLaZM2daPHbsWK9NP3O///57Su+dDGbaAAAAAAAAxBAPbQAAAAAAAGKIhzYAAAAAAAAxVKxr2miJb+eca9y4scXTp0/32g444ACLo/KKDz30UIt//PFHixOV6CuM5rEtX77cYi1nHVqwYIG3na9r2ihdvyJVN954o8X169dP2E9zSQvbRnzddNNNFoefGcZRZowYMcJiLcmdKi1tun79eq+tZs2aFmvZ2a+//trrt/3222/zfuS6MJ9byzbPmTPH4gceeCBr+3TyySdn7b2wpUaNGnnbzZo1S9hX723ef//9jO1TrqhWrZq3ffrppyfse+GFF1qs942ZpuvYjB49OmG/cE2bcD1IOHfDDTdYrCXckxWu03b88cdbHJYN1/VvMrkGRq6KWmemSZMmFmup59C4ceMs1t+V8+bN8/rVqFHDYl3L1Ln0rAOILenzgCuvvNLicIxVqFCh0L//6aefvO3PP//c4h9++MFr098gurZiixYtvH56TWjfvr3XNmnSJIu1bHi6MdMGAAAAAAAghnhoAwAAAAAAEEPFmh710UcfRW6rsFTbX8Jyo02bNrVYpzkdcsghSe/Xr7/+avGsWbMsDlO2dKqUTk3HtjnppJMs1tKZO+64o9dv2bJlFt96661e2y+//JKhvcO2qlWrlrfdvHlzi3W8OUdpxHQ5+uijve399tvPYp3em+xU33D6p05P1tKZzjnXtm1bi6PKEV9++eUW9+3bN6n9yDd33HGHt61TxHUqfpiilm763Rd+tpgunl1RKTuhMI0A0R555BFv++yzz7ZY7y+dc+7111/Pyj6FWrVqZfHuu+/utb3wwgsWv/zyy9napRJDU3edc+6CCy4otN/kyZO97aVLl1rcrl27hK9fsWJFizX1yjnnBg0aZPGSJUu2vrN5Lrz/Hzx4sMWaDuWcnx4clTKowpQoFS5/gfR7+umnvW1Na4sq363PDb777juLb7vtNq+f/q4PHX744RbrfeiAAQO8fvp8Qa8BzjnXp08fi4cOHWpxulNlmWkDAAAAAAAQQzy0AQAAAAAAiKFiTY9Kh9WrV3vbY8aMKbRfVOpVFJ16HKZi6VSsV199NaXXx5Y0XSacEqn0mH/66acZ3SekT5hOobJZdSPXaRraK6+84rVFTTdVWs1Lp3zec889Xr+odER9jUsuucTiqlWrev169uxp8c477+y1PfnkkxZv2rRpa7udUzp16mRxWLFg9uzZFmez0pqmuYXpUJ988onFP//8c5b2KH8dddRRCdvCqjRR6YnYUkFBgbetn/VFixZ5bZmsAFSmTBlvW6f+X3HFFRaH+9utW7eM7VMu0HQH55wrX768xVptJrxn0e+nf/7znxaHKRl169a1eI899vDa3nrrLYtPOOEEi1etWpXMrueFcuXKWRwugaDLKKxYscJre/jhhy1mqYT4CO/rtGrTRRdd5LWVKlXKYv1dEKbOP/TQQxanupxClSpVLNYqpnfffbfXT5dpCVMrs4WZNgAAAAAAADHEQxsAAAAAAIAY4qENAAAAAABADJX4NW0yoVq1ahY/9dRTFm+3nf+MS8tRk4eaujfffNPbPvbYYwvt99JLL3nbYflblAyNGjVK2KbrmmDb7LDD35f3ZNewCdeG6ty5s8Vh3niydE2bBx980OJevXp5/cqWLWtx+Dl4++23LZ4zZ05K+1FSnXHGGRbrMXLO/37KNF0jqWvXrhb/+eefXr/77rvP4nxbfyhbtESpxqEwx//bb7/N1C7lnRNPPNHb1nLqupZTuAZDsnQdldatW3tthx56aKF/88Ybb6T0Xvlqp5128rZ1TaBHH3004d9p+eDnn3/eYr1WO+dcnTp1Er6GrrWSyfWQSrJTTjnF4ltuucVr0zLcWvbeOefWrFmT0f1CasLr2I033mixrmHjnHM//fSTxbq27Ndff53Se+taNfvss4/Xpr8tR4wYYXG4jq0K93fgwIEWZ3ItP2baAAAAAAAAxBAPbQAAAAAAAGKI9KhCXHnllRZrWdqwvPjMmTOztk+5Zs8997Q4nN6tU1Y1JUOn3Tvn3Pr16zO0d0g3nc59wQUXeG0TJ060eNSoUVnbJ/yPlooOS8SmmhKViKY5aYqNc84dcsghaX2vkqpixYredqJUCOdST71IhZZr13S76dOne/3GjBmTtX3KV8mOlWx+PnJR7969ve02bdpYXL16da9NS6/r1PmOHTum9N76GmEpbzV37lyLw5LTiKblukOa/ham8CfSvHnzpN973LhxFnMvW7io1E+9b1y4cGE2dgfbSFOUnNsytVr98ccfFrds2dLiTp06ef3233//Qv9+48aN3vYBBxxQaOycf5+7++67J9wntXTpUm87W2nhzLQBAAAAAACIIR7aAAAAAAAAxBDpUc65I444wtsOVyn/i65k7pxzU6ZMydQu5byhQ4daXKVKlYT9Xn75ZYvzrWpMLmnXrp3FlStX9tpGjhxpsVZlQPqEle+UTj3NNJ3yH+5T1D7efffdFp9zzjlp3684CSua7LXXXhYPGTIk27tj6tatW+h/53sw+6LSMNJRuQj/M2HCBG+7cePGFjdt2tRrO/744y3WqijLly/3+r344otJvbdWI5k0aVLCfl9++aXF3CMVTXg91VQ2TUEMUzC0Auapp55qcVhtRsdi2HbxxRdbrOd62rRpyex6XghTYZSOt7vuustre+uttyymYl58fPzxx962plLrbwTnnKtRo4bFjz/+uMVRqaKabhWmYkVJlBK1efNmb3v48OEWX3PNNV7b4sWLk36/bcFMGwAAAAAAgBjioQ0AAAAAAEAM8dAGAAAAAAAghljTxjnXvn17b7t06dIWf/TRRxZ/9dVXWdunXKT5wgcffHDCfp988onFYa4qSqYmTZpYHOakvvHGG9nenbxw2WWXWRzm5haXDh06WHzQQQd5bbqP4f7qmja5bt26dd625uTrmhrO+etDrVq1Kq37Ua1aNW870foCY8eOTev7onBHHnmkxV26dEnYb82aNRZTCje9Vq9ebXFY2l63b7755m1+rzp16lisa4E5518Tbrjhhm1+r3w1evRob1vHjq5bE64zk2hdjfD1rrzySovfffddr23fffe1WNfH0O/tfFe1alWLw3sCXfvt3//+t9d2xx13WNyvXz+Ltcy6c/66KbNnz7Z46tSpCfepQYMG3rb+LuR6Gy0sw63rQe26665em64tq+vOrly50uu3YMECi/Uzob85nHOuRYsWRd7f/v37e9u33XabxbpeVTYx0wYAAAAAACCGeGgDAAAAAAAQQ3mbHlWmTBmLtXScc879/vvvFmt6zqZNmzK/YzkkLOWtU8s0BS2kU3/Xr1+f9v1Cduyxxx4Wt2rVyuKZM2d6/bSMHtJHU5GySac0O+fcgQceaLFeA6KEZXLz6dobTiHWMr6nn3661/bee+9Z3KtXryK/V8OGDb1tTcmoVauW15YoJSAuqXe5Tr9Pt9su8f9vGzVqVDZ2BxmmKR/h2NP0q/BaieSFKaVnnnmmxZq2XbFixYSv8cQTT1gcpsX9+uuvFg8bNsxr0/SP4447zuK6det6/fK5jPvDDz9s8fXXX5/03+n18Yorrig0Thcdf7q0Q+fOndP+XrksTDfS8ZGKl156yduOSo/SlHT9nL3wwgtePy0pXlyYaQMAAAAAABBDPLQBAAAAAACIIR7aAAAAAAAAxFDermlz4403WhyWnh05cqTFX375Zdb2Kdf861//8rYPOeSQQvu9+eab3jZlvnPD+eefb7GWD37//feLYW+QLbfffru3rWVPo8ybN8/i8847z2vTso75Rq+HYenfE0880eIhQ4YU+bVXrFjhbevaGbvttltSrxHmfSMzEpVcD9cCePrpp7OwN0i3M844w9s+99xzLdY1F5zbsuwt0kNLdut469Kli9dPx5yuPaRr2ITuvfdeb/uAAw6wuGPHjoW+nnNbfhfmE13X5NVXX/XaBg8ebPEOO/g/ZffZZx+Lo9b/Sgddw08/M1p23Dnn7rvvvozuB5y76aabLC7KmkKXXXaZxancR2UTM20AAAAAAABiiIc2AAAAAAAAMZQ36VE6jdw55+68806L165d67X16NEjK/uU65It0XfVVVd525T5zg01a9Ys9L+vXr06y3uCTBsxYoTF++23X0qvMW3aNIvHjh27zfuUK2bMmGGxlqR1zrmmTZtaXK9evSK/tpa1Db344ovedteuXQvtF5YoR3rsvffe3naYovGXhQsXetvjx4/P2D4hc0444YSEbe+++663/c0332R6d/KepkppnKrwOqnpPpoe1aZNG69f5cqVLQ5LlOc6LbEcXtfq16+f8O+OOeYYi0uXLm3x3Xff7fVLtGRDqjR9uVmzZml9bRTuoosuslhT0sKUOTV16lRve9iwYenfsQxhpg0AAAAAAEAM8dAGAAAAAAAghnI6PapKlSoWP/74417b9ttvb7FO7XfOuXHjxmV2x+DR6Z/OObdp06Yiv8aaNWsSvoZOj6xYsWLC19h111297WTTu3QK58033+y1/fLLL0m9Ri466aSTCv3v77zzTpb3JD/pVN2oCgpR0/L79+9vcfXq1RP209ffvHlzsrvo6dChQ0p/l8++/fbbQuN0mDt3blL9GjZs6G1PmTIlrfuRrw4//HBvO9EYDqsvomQKr8MbNmyw+JFHHsn27iDDXnvtNYs1Peqss87y+unyASzdkJyPPvqo0P+u6cTO+elRf/zxh8XPP/+81++ZZ56x+Nprr/XaEqWtIjNatGjhbeu1sVy5cgn/Tpfd0GpRzjn322+/pWnvMo+ZNgAAAAAAADHEQxsAAAAAAIAY4qENAAAAAABADOXcmja6Vs3IkSMtrl27ttdvzpw5Fmv5b2Tf5MmTt/k1Xn/9dW978eLFFu++++4Wh/nC6bZkyRJv+/7778/o+8XJkUce6W3vsccexbQncM65vn37WtyzZ8+E/bScbNR6NMmuVZNsv379+iXVD8VD10QqbPsvrGGTGbomX2jFihUW9+7dOxu7gwzQtRX0PsU555YtW2YxJb5zj35P6vfzySef7PW76667LH7llVe8tlmzZmVo73LThx9+6G3r/bmWiL744ou9fvXq1bO4devWSb3XwoULU9hDbE249mH58uUL7adrgjnnrxv1xRdfpH/HsoSZNgAAAAAAADHEQxsAAAAAAIAYyrn0qLp161rcrFmzhP20nLOmSiF9wlLq4bTPdDrjjDNS+jst8xeV1vH2229bPH78+IT9Pv/885T2Ixeceuqp3ramKk6cONHizz77LGv7lM+GDRtm8Y033ui1Va1aNWPvu3z5cm97+vTpFl9yySUWawoj4qegoCByG5l13HHHJWxbsGCBxWvWrMnG7iADND0qHF/vvfdewr/TlIBKlSpZrJ8LlBzffvutxf/+97+9toceesjiBx54wGs755xzLN64cWNmdi6H6L2Ic37Z9TPPPDPh37Vp0yZh259//mmxjtlbbrkllV1EIfR6d9NNNyX1N4MGDfK2P/nkk3TuUrFhpg0AAAAAAEAM8dAGAAAAAAAghnhoAwAAAAAAEEMlfk2bmjVretthSbe/hGs6aJlbZMZpp53mbWsuYunSpZN6jQYNGlhclHLdAwYMsHjevHkJ+w0dOtTiGTNmJP36+J+yZcta3L59+4T93njjDYs1BxiZM3/+fIs7d+7stZ1yyikWd+/ePa3vG5a579OnT1pfH9mx8847J2xj/YTM0O9FXZ8v9Ouvv1q8adOmjO4Tiod+T3bt2tVru+666yyeOnWqxeedd17mdwwZ9dJLL3nbl156qcXhPXWPHj0snjx5cmZ3LAeE31vXXnutxeXKlbO4efPmXr9q1apZHP6eGDhwoMV33333tu8knHP++Zg2bZrFUb8ddQzouc0lzLQBAAAAAACIIR7aAAAAAAAAxFCJT4/SErLOOVejRo1C+3366afeNuVLs69nz57b9PddunRJ054gXXRq/urVq702LZPeu3fvrO0TthSWWddtTSkNr6cdOnSwWM9n//79vX6lSpWyWKeyouS64IILvO2ff/7Z4nvvvTfLe5MfNm/ebPH48eO9toYNG1o8e/bsrO0TisdFF11k8YUXXui1PffccxYzFnPL8uXLve127dpZHKbm3HzzzRaHKXTYuqVLl1qs9zpaSt055w499FCL77nnHq9t2bJlGdq7/Na2bVuL9957b4ujfrtr2qimEOcSZtoAAAAAAADEEA9tAAAAAAAAYqhUUdKESpUqFYucoiOPPNLiESNGeG264rRq0aKFtx1OPY67goKCUlvvtXVxOYd5akJBQUHzrXfbOs5j8WEs5gTG4la888473navXr0sHjNmTLZ3p1C5PBarV6/ubd93330WT5gwweIcqM6Wt2NR72W1EpBzfgpr3759vTZNRf79998ztHdFk8tjMS7C6riHHXaYxS1btrR4G1KU83Ys5pJcGIuTJk2yuFGjRgn7PfTQQxZrumAOKHQsMtMGAAAAAAAghnhoAwAAAAAAEEM8tAEAAAAAAIihElnyu1WrVhYnWsPGOefmzJlj8fr16zO6TwAA5AotgYrsW7RokbfdrVu3YtoTZMrYsWMt1hK3QGE6derkbeu6H/Xq1bN4G9a0AWKhcuXKFpcq9fcSPWGJ9cceeyxbuxQLzLQBAAAAAACIIR7aAAAAAAAAxFCJTI+KotMFjznmGItXrVpVHLsDAAAAAClbu3att127du1i2hMgs3r16lVofO+993r9Fi9enLV9igNm2gAAAAAAAMQQD20AAAAAAABiiIc2AAAAAAAAMVSqoKAg+c6lSiXfGWlVUFBQauu9to5zWKwmFBQUNE/HC3Eeiw9jMScwFnMAYzEnMBZzAGMxJzAWcwBjMScUOhaZaQMAAAAAABBDPLQBAAAAAACIoaKW/F7hnJufiR1BpJppfC3OYfHhPJZ8nMPcwHks+TiHuYHzWPJxDnMD57Hk4xzmhkLPY5HWtAEAAAAAAEB2kB4FAAAAAAAQQzy0AQAAAAAAiCEe2gAAAAAAAMQQD20AAAAAAABiiIc2AAAAAAAAMcRDGwAAAAAAgBjioQ0AAAAAAEAM8dAGAAAAAAAghnhoAwAAAAAAEEM8tAEAAAAAAIghHtoAAAAAAADEEA9tAAAAAAAAYoiHNgAAAAAAADHEQxsAAAAAAIAY4qENAAAAAABADPHQBgAAAAAAIIZ4aAMAAAAAABBDPLQBAAAAAACIIR7aAAAAAAAAxBAPbQAAAAAAAGKIhzYAAAAAAAAxxEMbAAAAAACAGNqhKJ1LlSpVkKkdQbSCgoJS6XgdzmGxWlFQUFA1HS/EeSw+jMWcwFjMAYzFnMBYzAGMxZzAWMwBjMWcUOhYZKYNkD3zi3sHADjnGItAXDAWgXhgLALxUOhY5KENAAAAAABADPHQBgAAAAAAIIZ4aAMAAAAAABBDPLQBAAAAAACIoSJVjwKA4lKq1N8L4hcUsKg9AAAAgNzHTBsAAAAAAIAY4qENAAAAAABADJEehW2iKSvhdvny5b22bt26WXzqqadaXKZMGa9f6dKlLf7mm28srlSpktdvl112sXjEiBFe25AhQyxevXq1xb///nsh/wqkKpWUpfAzE9Wm29tt9/cz5j///NPrt3nz5qTeG0Wjx1zPxQ47+F8devw3bdqU+R3DFqLGFemEAABkDt/ByDRm2gAAAAAAAMQQD20AAAAAAABiiIc2AAAAAAAAMcSaNkirsmXLWty6dWuv7eyzz7a4QYMGFusaNs75eaFNmjRJ+F66dkadOnW8th9++MHiUaNGFfo3zqWWZ0reatHpMQvPd7Vq1Sxu3ry511ajRg2Lv/vuO4snTZrk9VuzZo3F4Xo3iLb99ttbrOtEOedcs2bNLL7hhhssrlWrltdP157q3bu31zZx4kSLOTfppeduxx13LPS/O+fcH3/8UWgcno+o61eiNY122mknr5+uG6bv5RxrT2VD1HiuUKGCxRs2bLB43bp1Xr/wvGHrdP0vHR/hGOMamFuK834w0XtzH5o5UfevDRs2tLh69epev++//97iOXPmeG16/6q/UcLvS64dYKYNAAAAAABADPHQBgAAAAAAIIZIj8I20SnBzjm38847WxxOD9Rphb/99pvFYcqSbutU73Cq4Pr16y3+4osvvDadiqjvxbTRbRc1HTgRPe7hFE+dir/33nt7bY0aNbJYp/BPmzYt6f3jnEfT4xNO/W3cuLHF9evXt7hcuXJev3333ddivQY4F12qHUUTXm933XVXi/fff3+LNVXKOefmz59v8ZIlSyz+9ddfU3pvTYPVfXDOubVr11r8yy+/eG16bY8al/k0ZvV6pXF4DBIdk/B6p+OvY8eOXttpp51msaY0PvPMM16/ZcuWbfV985Eea001c865ww47zOIqVapYPGHCBK/fggULLNZ7k/D+JpXjHl6/9TpQuXJlr03HoqZn6D4Vtl/5ItG4dG7L6/Bf0nEOo4Tvq+dbr8kbN270+hXlOp+v9LdGmTJlvDZNe7rwwgu9tk6dOllcvnx5i8PPjJ4D/Y50zr8mTJ8+3eJ+/fp5/caNG5f4H4C8wEwbAAAAAACAGOKhDQAAAAAAQAwVa3pU1JTDcHq3tmk6RZhak69TOeNCqyaEFWZ+/PFHi0ePHm3xhx9+6PX79ttvLdbzHlaS0mmKu+22m9emU0WzidScrQuPiY712rVre22aXjFz5kyLV69e7fXTcc8xT13dunW97auvvtpireQVlZLxz3/+02vTql86LZjzVHRhpSatwteuXTuLw+nxOl4WLVpkcbIpOM75Y0wrRIWfBU2dC6sQJZselcvCyl56/NJxHTvggAMsfvDBB702TWWrWrWqxf3790/pvXJdeK70+hhWydPUXk0900qWzvmpiipRyo1z0Wk3+ncVK1b0+mk1Rk3fcs65efPmWawVNhcvXpxwP/KJHldNe3HOuQMPPNBivSbrvatzftpZOn6bhNdarQ6n+6j34c751+tc/I0UlbKvbeF41vGi1W6vvPJKr1/Tpk0tDn9b6LGOGsOachX+vtXX1H3Ua7lzzv2///f/LCbVPHXhedJtPf5xfL7ATBsAAAAAAIAY4qENAAAAAABADPHQBgAAAAAAIIYysqZNmF+oZek0/1PLIjrn5+/VqVPHa9M8ec25DUv/rlixwmLN+Qvzw3UfwxxF/TvNGdXYOT9fdcOGDV7bypUrLdb1BcI8xJKe1x+WmNTztnDhQq9txIgRFmtuZlgWVul5CktWap6p5uc751z37t0tvvjiiy1ORx5oeM40HzLMlQzXdCiJUinxHSU8RrpWSr169by2qVOnWjxjxgyLw1zTkj6OipOOnWeffdZr0zWGovK1NZ++c+fOXtv69estfuyxxyxeunSp1y8O+cJxpMc9vM6dfvrpFmt59rfeesvr99NPP1mcjnVl9O/C9cT23HNPi2fNmuW1zZ071+J8Ot/pvoaqcP0KXY+hevXqCfdDvwvDErT5fD3VNQ32228/r+2pp56yWO8/nPNL+up6MVOmTPH6aUltPc5Rn5Fkz4euWeSccyeeeKLFen1wzrmxY8darOsK5vO512utrlHUo0cPr1/btm0t1uvYgAEDvH59+vSxeNWqVV5bKte/8DOiv6f0Ozi8p87k9ScOwn+fjmE9RuFvCP3+vO666yyuVq1awvfS9YHC7fD3kNLrw7p167w2XQvpmWeesThc7zNf17HR8xv1fEHvPRo1auT1O+WUUyw+9thjvTb9Dv35558tfvLJJ71+r7/+usXhb379HOh5Svd9DjNtAAAAAAAAYoiHNgAAAAAAADGUkfSosKyaloPV6Zv169f3+h133HEWh+WddQq2Tj1atmyZ10+nOe2+++4WV6hQweunU0DD6W76mvp6YblVTY/SkonOOTd06FCLdUp4SZ3elmh6mpaxc86fCqYpUM45N3HiRIt1inCy73vOOed4bZpWE6ZuaGnOqLSOdNDPUi6kQ4XC6Yi6rf/2qGnV+jfhOGrVqpXFOmad86cjUi46ffQcDBkyxGItIe1c4rETVSpaU1mdc+68886z+JhjjrH41Vdf9fpp2WHSNf6mYyec8tumTRuL9Rh9//33Xr9USs9GTanX7/iGDRt6bVpaeNiwYV6bpo3kq/C46hhLZSp1eG9z9NFHF/razvn3H6+88orF4T1QPtPPdrNmzbw2vecIUyEWLFhg8cCBAy0OUyEyeS3bZ599vG1N4QpTZsaNG2expu7kU9piSH9naDr//vvv7/XTMazppvobxjn/eqfH2zm/FHyyKavhudHvSR3ber0v7O9yQdQ9pZby1t+fYXrUXnvtZbHeuy9atMjrN3nyZIvHjBnjtWkKsJ7HmjVrev106Y6wNLyeR00nD5cByGVRaU/6HdeyZUuv39VXX22xPjcIlzPR63pYcl3Hh/7dLbfc4vU7/vjjLX777be9tnfeecdiPdekRwEAAAAAAOQBHtoAAAAAAADEEA9tAAAAAAAAYigja9qE67ZoyWvNWwvzBufPn29xWLpQ83E1/zZ8Ly17pyX7wrKYWtYrLMmoOcjNmze3WPOZnfPXZQlLm2qOouZKltS1GRLtt5axc84/Dpqb6VxqefN6zLt16+a1aY5imK997733Wqyfv0woqec0WcmuwRD2S1TONCwR/I9//MPicJxqme+Suh5UHITnpkuXLhYfeeSRFket/6TnOhzLeh0Ix4PmlO+7774W33rrrV6/gw8+2OJrrrnGa1u+fHnC/cp1mn/doUMHr61KlSoW6/dYmDOfKDc+HG+JSkKHypYta7HmeTvnXL169SzWssLO5ebaCsVBx6keb+e2LAuv9N5Jc/I5L3/T+wr9nIfC+wpd60LXt8n0/YGu3xHeI+n+Dx482GsbPXq0xeF9XL4I17bQY3TAAQdYHHVvo78XwrVkTjvtNIt1/THnnHv22Wct1vUfw+9Wfa9wnOp9r34e82E86zgNr3m6HoquL6Kxc8498sgjFn/00UcWh8fvu+++szgs3a6/7/TvwnspbYtaEzCf6LjS3+7OOXfmmWdafMEFF1gcroWrY1iPefibUM/bkiVLvLalS5darOsc6TXAOefatm1rcbhmka7Zqs8XwvuobT3XzLQBAAAAAACIIR7aAAAAAAAAxFBG0qPC6T86fUyn8C1evNjrpyWzNHbOn3aorx+WXdQpVlqmT6crOeen8YTTqKpXr25xjx49LNYy0s75/5bVq1d7bToVK9fKQOvxD0t369TOVNNZdIpb9+7dLdY0C+f8qY5hOsWHH35Y6P6i6FKdyqlTH/Wcatqic366VDi9eNmyZUnvZzL7ofLpcxFOPb3rrrssDq+hSq+bEyZMsPirr77y+mn6Ujil9JBDDrFYS7qHKbBayvH888/32nr16mVxvqXJ6XTdsCS7pqBquks4DVzpeNAp5s5Fp/Lq3+l+tGrVKuE+hdPM8+3c/SUqxSHRdSgqJUOvp40bN/b66XgOj7emYUR9RpIVVRa+pF5fdUyEKfH67w3TWDQlMZXPeVQ6Rdim35n9+vWzWFNdnXPu3XfftXjo0KFeW6bTxkuCsJS3fgfpuY5Kl3nooYcsDs/75ZdfbnE4TleuXGmxlpQuyjICOsZK6nhLVjgG6tSpY3FY6l7HYtT3kf720/So8LqWbEn2qOs8tjyueg943XXXeW0XXXSRxXq9i/pOe/755y3+4IMPvH56LxueQy3z3bVrV4v/7//+z+un361hmlbDhg0tnjRpUsL30mOQyphlpg0AAAAAAEAM8dAGAAAAAAAghjKSHhXSKUA6/VpXXXfOT42Iml6aqCqNc/40fa1GFU5Vi0pZ0rQtnXYXrjSv0xg//vhjry1fVuNPxxTA8BxqFZn27dtbHB5TTZl44403vLZ8nYafCeEUvkSr4EdNA9SprR07dvT67bnnnhZrtSjn/CnEyU4lDKfR5vq04USijnmYovaXsOLbbbfdZrGOsQoVKnj99FwPGTLEa9N00zvvvNNirRrmnP85CNOjevfubXGuj+0wZUkrboXmzZtn8ZtvvmlxVAWSqGn/UcdWK01df/31FofV4HTKuabUFfZ++SjZ6fUhPW877bSTxZdcconXT+9Tws/BE088YXGyadtRKVC5KKoioh6zMF1bK4toqmKYOq/0NcIUVh1Hem10zrnnnnvO4kMPPdTiMOVJ++VzBT6l34thuqlee/VY9u/f3+t3//33F9rvhhtu8Po1adKk0PcN3ytf71GKonbt2t726aefbrF+Dzrnj+Goqk1Kv/uiqmiGOHfJ0zQk55xr3bq1xRdeeKHXpqn0em7CipTnnXeexWFVKJXsd6tWVg2v8Yn+xjl/qZCo+xyqRwEAAAAAAOQgHtoAAAAAAADEEA9tAAAAAAAAYigra9oozefSMmphWyqv55yf+6ZxsvlszvlrCGhZrzBPTctKf/bZZ14bufvRNGdUS6U559wDDzxgseZ5jxs3zus3YMAAi6OOdy6WJc2mVNdgUHoejzrqKK9NPwvhOEq2/KWeY873/1SqVMni++67z2vTY65j55NPPvH66ZoIup6Dlk8MhWtlaN9BgwZZHH4OypYtW+i+O+eXhkxHGfg4K1OmjLfdrFkzi8PP78CBAy1eunRpwn4qUb7/1ujaRC1atEj4Xu+8847Fc+bMSfr180Wy19PwOqZjVksVh2WL9e/CdUzGjx9vcbLnPt+up3pfGq6xpus6heszaLntHj16WKzHPHz9qlWrWhxeU/W9dG0x55yrWbNmofv+zTffeNtff/21xbl4rlKhx7VWrVpe2w8//GCxruH28MMPe/10vSF9DS3x7Zx/36NrXjjnf08mu+5KvtF1f6699lqvTdfE099izvlrSiV7vdU4XNstn0qrp5t+f1SuXNlru+aaayyuVq1awtfQtRafeeYZr02/4/S9wrUBo9p0faTGjRsn3A8V/jZZsWKFxZn8jDDTBgAAAAAAIIZ4aAMAAAAAABBDWU+PUpmeZpbs65cuXdrbPvfccy3WKfth6cY+ffpYHJZaxJb0OLds2dLiRx991OtXt25di3XKfzgtbu3atQnfK2oqnEq27Cn+luy40nOg5RrDctOLFi2yeOTIkV5bOs5PvkxtDctUailYLZ8Y0mn53bp189oSpacVJWVO0wEmTpxocTgudf/D8rdaYjwX06N0rFSpUsVr03K/YbrLV199ZXG6r2Xh5+nwww+3WEtO6zXaOf96HqZAY0uJxlL43/X7U8uS6j2Kc/7Ufi0D71x0+elEwvSoRGlVuXJt1X/ftGnTvLaZM2da3Lx5c69NUzh1un3Hjh2Tev3wuqblovfZZx+vTc+Jpg7ccsstXj/G35b0uqbpUM756aaaHhWmNmn67uOPP27xbrvtlvB9w9dYuXKlxTq2w/GWK+MqFXvssYfFZ511ltem6YnhWIy6508kn49zJul4q1GjhtemKVFR3zOrVq1K+Prt27e3WFPuw2U3dKzr/YtzzrVr187iHXfcMeF76Xfr7Nmzvbbvvvsu4d+lEzNtAAAAAAAAYoiHNgAAAAAAADFUrNWj0iGcwp1K6oZWiHLOueOPP95inXIerlCuU2WpFrWlcIqipkQNHjzYYp3+75w/pVeP+aeffur109SN8LxrhYCdd97Z4nCKKtInnN6oUxCPOOIIi8N0xFGjRlmsqVLOpTaew7GYL9Ne9TPvnJ/qFLbp2Ln33nst1inbUVI9ppr6E34O9Bz++uuvXlu4nWv03x5WztLUiPA4pDslSvdD0z2cc+6iiy4q9G+GDx/ubS9evDit+5SvwuupVt5o27Ztwn6aNqzV35zbsipKMu+db9dT/feG1c+0+l1Y7UnHbZiypvT8DBkyxOIwlalOnToJX0PHvY6/b7/9NuHf5Kuo6mdR1ypN0z/66KO9tksvvdTiBg0aWByOlTVr1lg8d+5cr22//faz+LjjjrNYq+8551eqyjdaySf8PtLfF2EFPU3H17SV8NoV/n5M9N+jqnvp50v/LrzW5utvxKjvEr3mhan4eq+zbt06i2+99Vavn/5+D9OelF4zw2uCbuv5De+3Jk2aZHHnzp29NqpHAQAAAAAA5DEe2gAAAAAAAMQQD20AAAAAAABiqFhLfhdForzBMPcw2Rx/Xefkqquu8tq0DNmPP/5occ+ePb1+ub7OwrbSknzOOffss89aXL169YR/N3/+fIt79eplcVjiO9myw1qOPV/zSrMhag0GLT8driv09ddfW5yoxPTW3gtbHhO9NoZjZcOGDRa//PLLFic7PqKOf9T6UieffLLF4boP+nfkg/9N1/4Jy4G3adPG4tdff91iveaF9HMRriuk5+Tss8/22nQM6/n45ptvvH7JrpuCaOEY0zVO9tprL4vD8bZw4UKL582bl9J7Jcrxzwf67w3vObQM9Pvvv++16VjS8ReOB13HQWO9J3XOL1nbtGlTr03X09F7pGS/P0P5dL71fIT38K1bt7b4mGOOsXjvvff2+um51jVnwnUv9TNy4IEHem2HH364xQcddJDF4feiliFP9vdNSS4brvteq1Yti6PWkilXrpzX9sADD1h8//33F/o3zjl32GGHWay/SVavXu31++CDDywO17LSa7GO52nTpnn9ospW5zIdb5MnT/baLrvsMou7dOnitemx1PLaLVq08Prpmn96rxmu46ilvMPPko4rXefqySef9Po98cQTFofXjmzdozLTBgAAAAAAIIZ4aAMAAAAAABBDsU2PCqexJZr2FE4HTTQNMEyj+sc//mFxx44dvTZ9TZ0ONX369KTeC/+z7777ets6jVDPRzjdsGvXrhZrelpRjndUib5tVZTUnJLyGSlKukuivwtfQ0u5a1nMcFqhjqtkUyvC8az7WFKOebqF5Q512mg4dVPLE4alZhPR8xu+V9Tx12nhl19+ucVasjPcRy3T6Vx+TS0Op2brNTCcGnzmmWdarCVqP/30U6+fpmvoNPBwWrkKS3yXKVPGYh3D4bnJ1/GnUr2eqvAad+SRR1pcoUIFi8Ox/dlnn1msaZBRotIp8vl8ht9HOo6iUhBTEd7LaipkmFL89ttvWzxz5kyLk/2udi5x+myun289DuFx1XQcTYnS1Arn/Ovfgw8+aHHfvn0Tvu+5557rbXfo0MFiLWetKSPO+ed65cqVCV9fleRzqOdHr3Ph+NDfgeG/V0uFDxo0qNDXC7f1fcNxf+GFF1ocft/pONKlHW6//XavXz7dwyg9N+H3kS6NEKZZJ0o3eu6557ztf//73xbrGLv66qu9flWrVi10n5xzbsmSJRafc845Fn/xxRdevzikfjPTBgAAAAAAIIZ4aAMAAAAAABBDPLQBAAAAAACIoRKzpo3moGkearJltrT8sHPO3XnnnRaHJfZGjx5t8eDBgy2OQz5b3GkpxJtvvtlr0zUR9Fi++eabXj9dzyLVEsTpyOnVf4uuvxGuNaC5tWHOZpw/M1Hr0SR7/LRf+BqNGjWyWNe3mTVrltcv2XWL9LiH75XPJaH/onnxzjlXs2ZNi8Pyh+XLly/y6+sYCHP89dyEZam1ZGnYpnSsjBo1ymsL1x7INfq5X758udemZWP1Guqcv+5Cy5YtLdby0M755YT1WA4fPjxhv1122cVr0zEXVVozn8oHq3RcT1U4xjp16lRoW7i2it6zRJUIDr/HVD6dt7gI1wDUNU/WrFnjtT366KMW6/oqRVnTRsetfk5KcrnoogqvtVraWcdYeF/3n//8x2JdlyMcb/qdGa6PcdZZZ1ms19pwvThdZ2fdunVem65Hl4vnaerUqRYvWLDAa9ttt90sXrt2rdem94O6lknFihW9fomugeEY0PuW8B5GPyd6rsLxPGXKlELfK5+En1HdDtcsStbChQstfvXVVy0+++yzvX76eQnfS78z//vf/1ocx99vzLQBAAAAAACIIR7aAAAAAAAAxFBs06PCdIdU0h80veXiiy/22rQ8apiucffdd1scTrtDtAMOOMBiLVHqnD8VUafo6/RS5/xzHTXlXKeeJjuNLWqKcI0aNby21q1bW6zTaMOyizpFNfwsheV74ypq2mKydLw559xRRx1lsZYWDqeJhlN+VdT5hy8s36zTdsNjp9OxK1WqZHF4LvRzoK+/1157ef20jPT555/vtek1Ieocfv/99xY/88wzCfcjF+m/L0x3ee+99yz+/PPPvTZNj9I0X02Nc84fmzNmzLBYp5+Hr3Haaad5bbVr1y5036PSqHL9vCWSjn+3Tut3zp96r8d40aJFXr/Jkycn3A+uodsmHZ9tfQ0dOw8//LDXT9NdBwwY4LVpaeFUUpmdS/4+K9fGsN4rLl261Gt76qmnLH7yySctDo+BvkbUmNJjHN733HHHHRZfddVVFut9rXP+fWmYIqT3l7pPJfmc6b5rqsr111/v9dNU3vA86j3HCSecYPFJJ52U8DVU+HtTf7uE6cD63ar9WrRo4fV76623Er4+UqfHX9NGwxRx/VyNHz/ea3v66actTjVNK1uYaQMAAAAAABBDPLQBAAAAAACIodimR6WDTo+64oorvDadFvfaa695bTNnzrS4JE8zzIZwamirVq0sDit2Jao+olVPnHNu5MiRFus0wvr163v9GjZsaPGYMWO8tmXLllmsq/GH0yNPPvlkixs3buwS0Ypi06ZN89pWrVplcbhi/Msvv5zwNYtbuj/b1apV87bbtWtXaL8vv/zS245KbdN9ZCxG08+8c37aRL169bw2HRNa3enGG2/0+ulnu0uXLhafccYZXj+91obTh3W6t47nsHLH6aefbvHPP//s8lX4OdfpumFqpk6P13S4cCq+jjFNvwqnaet7aQpG+Bqa3hrHCgvFIR3XJ/2OPOigg7w2TaXRYz5ixAiv3y+//JLw9bmGFk1U6ktU9S0VlaK2xx57WNykSROvn143w7TIqKpgqdDrQK5/RqKuV6kc16jPgR7XMO1Cq0lpek+4rID+Vtl99929Nk3NX79+vcXhv7EknVPdV/2++/jjj71+USlqEydOLDQOU/jbtGljsVZmDI+X3tOE6Wv63nq+9dyEf0d6VOrCdOw33njDYv3NEZ4nrbLXs2dPr+2nn36yOO5jhZk2AAAAAAAAMcRDGwAAAAAAgBjioQ0AAAAAAEAM5dyaNpqX2Lt3b4vDXNB58+ZZPGjQIK8tzEVE8tasWWNxVD645huGpWXbtm1rcfny5S3Wc+ucnxeqa28456+JoaVTK1So4PXTfORwfzVHWP8uXDtkv/32s/jbb791+USPn6554pxfBlrzucP1NpLN7417rmlxC9c7efDBBy3u06eP16ZjSUtTanlp5/wxoWMxKq87pLnnuk5K165dvX6zZs2ymHNduKjSs5qzrXFhf5eIjkW9ljvnr2Oj18aFCxem9F7Yko4rXePJOX89Bj0XgwcP9volu8aQnidKgScn2dLYUWux6d/pGnjlypXz+ul35o8//pjiHieW6Hs3F8ZvpkuYJ3r98JhGvZf+ztB7yjlz5nj9dK0yjZ3zv8f1mp8r64zp8SzKOjD67588ebLF99xzj9dPS4PXrFnT4nCdIn3vqM/T2rVrLQ6/F1P5TGbis1sS6RqMb775ptemvxf1vIW/43v06GGxrpvqXPrXCMskZtoAAAAAAADEEA9tAAAAAAAAYqjEp0eF5WW1ZG3r1q0T/p2mDmi5PRRNOF3vo48+snjq1KleW4MGDSzW86ZpF85tmcKUiE5ZDNPfqlSpYnGYyqF0WpyWfXPOuRkzZlj84osvWhz+u95///2Er5HrdEp3+/btvTadqrhixQqL8+0YZUs4FrUUcJi2pylROhYrVark9Us2bULfW1M3nHPus88+s/i6666zWNOhnKMM5rZKx9RpvVZqeqNz/nldsmSJxWFqKlKn19MDDzzQa9Pzq9dTTfUuClKiti4cU+kYY5rmduGFF1pctmxZr9+GDRss1hRv5/xxGjW1PyolQ7dzIe1C/z1RJZr1eKUjLSLVFB7dX017irr/DVNW86lUe6r0GOk9vXPOTZs2zeJ99tkn4Wv88ssvFof3N7o9d+5ci5cvX+71089kmLoTNU5VPp1j/f1w6aWXWnz00Ucn7Kfn+umnn/b6/ec//7G4JB9HZtoAAAAAAADEEA9tAAAAAAAAYoiHNgAAAAAAADFUIte00Zy//fff32vTHGHNHR4/frzXb+jQoRazlkL66PpA7dq189o6dOhg8QUXXGCxlr10zrnKlStbHJUHunr1aosXLVrktWn5Yy2LqGsxOOfcu+++a/HXX3/ttem/ZePGjRaH5RSjynvmIh1/eq7C/HA9P1piT8uxp/q+oXw47kWlY+D888/32vr372+xrm+jY8U5P19Yj384BnSdoksuucRr++STTyzWMracs+IXjildXyxs0/UUtCxtWF4cRaPHWddmC8vO6nGeP39+of891fdlLGZOOI723HNPi4844giLw7VXdJ2TQw45xGv78ssvLV6/fr3FRbk30b65cP7136D3/rvttpvXT9cg0RLNYVuyxyTV3w96vnUdm/D1dKzr2irOJb4vzYXzmQl6/+Gccy+99JLFNWrUsDi89upaNbpup3P+2lP6vThz5kyvn56rqPPDufsfXVPviiuusDi8Tup4GT16tMXdu3f3+qVyXMPPQRx+6zHTBgAAAAAAIIZ4aAMAAAAAABBDJTI9SqeNnnrqqV6bTjPUqYQvvPCC10+ntCF9dKpaWApWy2ZrHCXZkpXJCv8mnE6sEk1/C18jn6cz6hh77bXXvLavvvqq0LYwzS1Z+XycU6HHS6ftOufcCSecYHH9+vUtPvLII71+tWvXtljLDI8aNcrrN336dIvDkpjYNplMY4lKjwrfS6d3a5lp/T4OX5Mxu3V6vPT4h2mkWkL2hx9+sDicwh11/Dk3xW/nnXe2WFPbwrQYPVeaUuWcn5a8bt06i8Nzqq+ZTyXe9TtI08eci05F0mOk9ylRaWfJCo+/pv7rtTUc97ofYVl4/TsVXhNYAqJwH3/8scUTJkyweJdddvH6aRpdeH70vOrnLh2fmXyiKY3OOXf66adbvMcee1gcfpZ1aYyLL744Yb9UFCU9KlupU8y0AQAAAAAAiCEe2gAAAAAAAMRQiUmP0hWjW7VqZfG5557r9dMphzqdWFM1nGO6YEkRh9W6Q3HZj+Ki/35NmRkwYEDCflpZKtWxl+/HPZ20isKUKVMKjREPmfzch6+t35mDBw/22jTdVasx6t84l19pGOmg18MZM2ZYfOutt3r9ypUrZ/Hs2bMtDivgUJkkXsJjvnDhQotfeeUVizt16uT102n/w4cP99q0CqZ+t3Lu/0dTU7TqnXN+StpOO+3ktelYTBSnKjz+mvak19Bwf7WiY5hCoq+hr8/vm8KFx0WvneF1FJmn9woVK1b02tq2bWuxjtOwWmK/fv0sDqsIpyJRxdSwLfwsRS21kU7MtAEAAAAAAIghHtoAAAAAAADEEA9tAAAAAAAAYihWa9po/piuYeOcc7Vq1bL40UcftbhevXpeP80zmzt3rsVaotS5/MrvBTJFxxulnoGSJfwe/OWXXyweO3as1zZx4kSLNX97w4YNXj/WU0idlif++uuvvTa9P8pW/jzST8dLjx49LB46dKjXT9crmTlzptcWruuAxMKxotc4XdvNOf96mOkxlmidv3CfdH91zc7wNfhNg5JGf+d369bNa0u0po1eF51zbsyYMRbreEiH8PXiMMaYaQMAAAAAABBDPLQBAAAAAACIoVilR6mw1FblypUt/vnnny0Op4nq1EItS7pu3bo07yEAALlDp/+G05D1exeZR5pZbtIxtnHjRos1/dA5UuCyITzG25r+EP5uSfb1osp1a0pUeE0Ot4GSRNOPhg0b5rUdf/zxFu+5554W9+3b1+v3zTffpHWf4v69y0wbAAAAAACAGOKhDQAAAAAAQAzx0AYAAAAAACCGShUlh7NUqVLFVu9qu+3+fr5UunRpi8MSeJqjmktlEQsKCkptvdfWFec5hJtQUFDQPB0vxHksPozFnMBYzAGMxZzAWMwBjMWcwFjMAYzFnFDoWGSmDQAAAAAAQAzx0AYAAAAAACCGilrye4Vzbn4mdmRrtAzXb7/9Vmicw2qm8bWK7RyC85gDOIe5gfNY8nEOcwPnseTjHOYGzmPJxznMDYWexyKtaQMAAAAAAIDsID0KAAAAAAAghnhoAwAAAAAAEEM8tAEAAAAAAIghHtoAAAAAAADEEA9tAAAAAAAAYoiHNgAAAAAAADHEQxsAAAAAAIAY4qENAAAAAABADPHQBgAAAAAAIIb+P+AdluwyllY3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x288 with 20 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "n = 10  # How many digits we will display\n",
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    # Display original\n",
    "    ax = plt.subplot(2, n, i + 1)\n",
    "    plt.imshow(x_test[i].reshape(28, 28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "    # Display reconstruction\n",
    "    ax = plt.subplot(2, n, i + 1 + n)\n",
    "    plt.imshow(roundtrip_imgs[i].reshape(28, 28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b263f25",
   "metadata": {},
   "source": [
    "What about the compression? Let's check the sizes of the arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "fb73c4dc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:49:40.665212Z",
     "iopub.status.busy": "2021-11-29T20:49:40.664326Z",
     "iopub.status.idle": "2021-11-29T20:49:41.288376Z",
     "shell.execute_reply": "2021-11-29T20:49:41.289220Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_test size (in MB): 29.91\n",
      "encoded_imgs size (in MB): 1.22\n",
      "Compression ratio: 1/25\n"
     ]
    }
   ],
   "source": [
    "encoded_imgs = autoencoder.transform(x_test)\n",
    "print(f\"x_test size (in MB): {x_test.nbytes/1024**2:.2f}\")\n",
    "print(f\"encoded_imgs size (in MB): {encoded_imgs.nbytes/1024**2:.2f}\")\n",
    "cr = round((encoded_imgs.nbytes/x_test.nbytes), 2)\n",
    "print(f\"Compression ratio: 1/{1/cr:.0f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb9c172a",
   "metadata": {},
   "source": [
    "## 6. Deep AutoEncoder"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da51e981",
   "metadata": {},
   "source": [
    "We can easily expand our model to be a deep autoencoder by adding some hidden layers. All we have to do is add a parameter `hidden_layer_sizes` and use it in `_keras_build_fn` to build hidden layers.\n",
    "For simplicity, we use a single `hidden_layer_sizes` parameter and mirror it across the encoding layers and decoding layers, but there is nothing forcing us to build symetrical models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3097ac1d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:49:41.293867Z",
     "iopub.status.busy": "2021-11-29T20:49:41.292603Z",
     "iopub.status.idle": "2021-11-29T20:49:41.303017Z",
     "shell.execute_reply": "2021-11-29T20:49:41.303880Z"
    }
   },
   "outputs": [],
   "source": [
    "from typing import List\n",
    "\n",
    "\n",
    "class DeepAutoEncoder(AutoEncoder):\n",
    "    \"\"\"A class that enables transform and fit_transform.\n",
    "    \"\"\"\n",
    "    \n",
    "    def _keras_build_fn(self, encoding_dim: int, hidden_layer_sizes: List[str], meta: Dict[str, Any]):\n",
    "        n_features_in = meta[\"n_features_in_\"]\n",
    "\n",
    "        encoder_input = keras.Input(shape=(n_features_in,))\n",
    "        x = encoder_input\n",
    "        for layer_size in hidden_layer_sizes:\n",
    "            x = keras.layers.Dense(layer_size, activation='relu')(x)\n",
    "        encoder_output = keras.layers.Dense(encoding_dim, activation='relu')(x)\n",
    "        encoder_model = keras.Model(encoder_input, encoder_output)\n",
    "\n",
    "        decoder_input = keras.Input(shape=(encoding_dim,))\n",
    "        x = decoder_input\n",
    "        for layer_size in reversed(hidden_layer_sizes):\n",
    "            x = keras.layers.Dense(layer_size, activation='relu')(x)\n",
    "        decoder_output = keras.layers.Dense(n_features_in, activation='sigmoid', name=\"decoder\")(x)\n",
    "        decoder_model = keras.Model(decoder_input, decoder_output)\n",
    "\n",
    "        autoencoder_input = keras.Input(shape=(n_features_in,))\n",
    "        encoded_img = encoder_model(autoencoder_input)\n",
    "        reconstructed_img = decoder_model(encoded_img)\n",
    "\n",
    "        autoencoder_model = keras.Model(autoencoder_input, reconstructed_img)\n",
    "\n",
    "        self.encoder_model_ = BaseWrapper(encoder_model, verbose=self.verbose)\n",
    "        self.decoder_model_ = BaseWrapper(decoder_model, verbose=self.verbose)\n",
    "\n",
    "        return autoencoder_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ddab7a22",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:49:41.308046Z",
     "iopub.status.busy": "2021-11-29T20:49:41.306825Z",
     "iopub.status.idle": "2021-11-29T20:50:35.232982Z",
     "shell.execute_reply": "2021-11-29T20:50:35.231752Z"
    }
   },
   "outputs": [],
   "source": [
    "deep = DeepAutoEncoder(\n",
    "    loss=\"binary_crossentropy\",\n",
    "    encoding_dim=32,\n",
    "    hidden_layer_sizes=[128],\n",
    "    random_state=0,\n",
    "    epochs=5,\n",
    "    verbose=False,\n",
    "    optimizer=\"adam\",\n",
    ")\n",
    "_ = deep.fit(X=x_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0adbb576",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-11-29T20:50:35.239122Z",
     "iopub.status.busy": "2021-11-29T20:50:35.235960Z",
     "iopub.status.idle": "2021-11-29T20:50:36.696673Z",
     "shell.execute_reply": "2021-11-29T20:50:36.697572Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1-MSE for training set (higher is better)\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AutoEncoder: 0.9899\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Deep AutoEncoder: 0.9902\n"
     ]
    }
   ],
   "source": [
    "print(\"1-MSE for training set (higher is better)\\n\")\n",
    "score = autoencoder.score(X=x_test)\n",
    "print(f\"AutoEncoder: {score:.4f}\")\n",
    "\n",
    "score = deep.score(X=x_test)\n",
    "print(f\"Deep AutoEncoder: {score:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf1f3736",
   "metadata": {},
   "source": [
    "Suprisingly, our score got worse. It's possible that that because of the extra trainable variables, our deep model trains slower than our simple model.\n",
    "\n",
    "Check out the [Keras tutorial](https://blog.keras.io/building-autoencoders-in-keras.html) to see the difference after 100 epochs of training, as well as more architectures and applications for AutoEncoders!"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}