Course_4.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZMi612g2a1wO"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/neworldemancer/DSF5/blob/master/Course_4.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Nd5wH914a1wU"
      },
      "source": [
        "# Introduction to machine learning & Data Analysis\n",
        "\n",
        "Basic introduction on how to perform typical machine learning tasks with Python.\n",
        "\n",
        "Prepared by Mykhailo Vladymyrov & Aris Marcolongo,\n",
        "Data Science Lab, University Of Bern, 2023\n",
        "\n",
        "This work is licensed under <a href=\"https://creativecommons.org/share-your-work/public-domain/cc0/\">CC0</a>.\n",
        "\n",
        "# Part 4."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LzEnLAMt3De2"
      },
      "source": [
        "# 0. Load libraries"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hVJn0ilgOS8F"
      },
      "outputs": [],
      "source": [
        "from matplotlib import  pyplot as plt\n",
        "import numpy as np\n",
        "from imageio import imread\n",
        "import pandas as pd\n",
        "from time import time as timer\n",
        "\n",
        "import tensorflow as tf\n",
        "import plotly.graph_objects as go\n",
        "\n",
        "%matplotlib inline\n",
        "from PIL import Image\n",
        "from matplotlib import animation\n",
        "from matplotlib import cm\n",
        "from IPython.display import HTML\n",
        "from sklearn.datasets import make_blobs\n",
        "\n",
        "plt.rcParams[\"animation.html\"] = \"jshtml\"  # for matplotlib 2.1 and above, uses JavaScript"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "joJSww74bWnx"
      },
      "outputs": [],
      "source": [
        "import os\n",
        "import tarfile"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7ZPG2uLJbSju"
      },
      "outputs": [],
      "source": [
        "if not os.path.exists('data'):\n",
        "    path = os.path.abspath('.')+'/colab_material.tgz'\n",
        "    tf.keras.utils.get_file(path, 'https://github.com/neworldemancer/DSF5/raw/master/colab_material.tgz')\n",
        "    tar = tarfile.open(path, \"r:gz\")\n",
        "    tar.extractall()\n",
        "    tar.close()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tWsA0mkObhN-"
      },
      "outputs": [],
      "source": [
        "from utils.routines import *"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pclZR6uFklf_"
      },
      "source": [
        "# 1. Datasets"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s_wxOrdWko8W"
      },
      "source": [
        "In this course we will use several synthetic and real-world datasets to illustrate the behavior of the models and exercise our skills."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8UQgU5I-lEll"
      },
      "source": [
        "## 1. Synthetic linear"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jGfWOWRjlWPa"
      },
      "outputs": [],
      "source": [
        "def get_linear(n_d=1, n_points=10, w=None, b=None, sigma=5):\n",
        "  x = np.random.uniform(0, 10, size=(n_points, n_d))\n",
        "\n",
        "  w = w or np.random.uniform(0.1, 10, n_d)\n",
        "  b = b or np.random.uniform(-10, 10)\n",
        "  y = np.dot(x, w) + b + np.random.normal(0, sigma, size=n_points)\n",
        "\n",
        "  print('true slopes: w =', w, ';  b =', b)\n",
        "\n",
        "  return x, y"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5RLYxGy_nBZG"
      },
      "outputs": [],
      "source": [
        "x, y = get_linear(n_d=1, sigma=1)\n",
        "plt.plot(x[:, 0], y, '*')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "10ODDOp4nX4S"
      },
      "outputs": [],
      "source": [
        "n_d = 2\n",
        "x, y = get_linear(n_d=n_d, n_points=100)\n",
        "\n",
        "fig = plt.figure(figsize=(8,8))\n",
        "ax = fig.add_subplot(111, projection='3d')\n",
        "ax.scatter(x[:,0], x[:,1], y, marker='x', color='b',s=10)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FJ5rjq7fIe8Q"
      },
      "source": [
        "## 2. House prices"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A-45usskInlD"
      },
      "source": [
        "Subset of the Ames Houses dataset: http://jse.amstat.org/v19n3/decock.pdf"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dVv2ID96IyN0"
      },
      "outputs": [],
      "source": [
        "def house_prices_dataset(return_df=False, return_df_xy=False, price_max=400000, area_max=40000):\n",
        "  path = 'data/AmesHousing.csv'\n",
        "\n",
        "  df = pd.read_csv(path, na_values=('NaN', ''), keep_default_na=False,  )\n",
        "\n",
        "  rename_dict = {k:k.replace(' ', '').replace('/', '') for k in df.keys()}\n",
        "  df.rename(columns=rename_dict, inplace=True)\n",
        "\n",
        "  useful_fields = ['LotArea',\n",
        "                  'Utilities', 'OverallQual', 'OverallCond',\n",
        "                  'YearBuilt', 'YearRemodAdd', 'ExterQual', 'ExterCond',\n",
        "                  'HeatingQC', 'CentralAir', 'Electrical',\n",
        "                  '1stFlrSF', '2ndFlrSF','GrLivArea',\n",
        "                  'FullBath', 'HalfBath',\n",
        "                  'BedroomAbvGr', 'KitchenAbvGr', 'KitchenQual', 'TotRmsAbvGrd',\n",
        "                  'Functional','PoolArea',\n",
        "                  'YrSold', 'MoSold'\n",
        "                  ]\n",
        "  target_field = 'SalePrice'\n",
        "\n",
        "  df.dropna(axis=0, subset=useful_fields+[target_field], inplace=True)\n",
        "\n",
        "  cleanup_nums = {'Street':      {'Grvl': 0, 'Pave': 1},\n",
        "                  'LotFrontage': {'NA':0},\n",
        "                  'Alley':       {'NA':0, 'Grvl': 1, 'Pave': 2},\n",
        "                  'LotShape':    {'IR3':0, 'IR2': 1, 'IR1': 2, 'Reg':3},\n",
        "                  'Utilities':   {'ELO':0, 'NoSeWa': 1, 'NoSewr': 2, 'AllPub': 3},\n",
        "                  'LandSlope':   {'Sev':0, 'Mod': 1, 'Gtl': 3},\n",
        "                  'ExterQual':   {'Po':0, 'Fa': 1, 'TA': 2, 'Gd': 3, 'Ex':4},\n",
        "                  'ExterCond':   {'Po':0, 'Fa': 1, 'TA': 2, 'Gd': 3, 'Ex':4},\n",
        "                  'BsmtQual':    {'NA':0, 'Po':1, 'Fa': 2, 'TA': 3, 'Gd': 4, 'Ex':5},\n",
        "                  'BsmtCond':    {'NA':0, 'Po':1, 'Fa': 2, 'TA': 3, 'Gd': 4, 'Ex':5},\n",
        "                  'BsmtExposure':{'NA':0, 'No':1, 'Mn': 2, 'Av': 3, 'Gd': 4},\n",
        "                  'BsmtFinType1':{'NA':0, 'Unf':1, 'LwQ': 2, 'Rec': 3, 'BLQ': 4, 'ALQ':5, 'GLQ':6},\n",
        "                  'BsmtFinType2':{'NA':0, 'Unf':1, 'LwQ': 2, 'Rec': 3, 'BLQ': 4, 'ALQ':5, 'GLQ':6},\n",
        "                  'HeatingQC':   {'Po':0, 'Fa': 1, 'TA': 2, 'Gd': 3, 'Ex':4},\n",
        "                  'CentralAir':  {'N':0, 'Y': 1},\n",
        "                  'Electrical':  {'':0, 'NA':0, 'Mix':1, 'FuseP':2, 'FuseF': 3, 'FuseA': 4, 'SBrkr': 5},\n",
        "                  'KitchenQual': {'Po':0, 'Fa': 1, 'TA': 2, 'Gd': 3, 'Ex':4},\n",
        "                  'Functional':  {'Sal':0, 'Sev':1, 'Maj2': 2, 'Maj1': 3, 'Mod': 4, 'Min2':5, 'Min1':6, 'Typ':7},\n",
        "                  'FireplaceQu': {'NA':0, 'Po':1, 'Fa': 2, 'TA': 3, 'Gd': 4, 'Ex':5},\n",
        "                  'PoolQC':      {'NA':0, 'Fa': 1, 'TA': 2, 'Gd': 3, 'Ex':4},\n",
        "                  'Fence':       {'NA':0, 'MnWw': 1, 'GdWo': 2, 'MnPrv': 3, 'GdPrv':4},\n",
        "                  }\n",
        "\n",
        "  df_X = df[useful_fields].copy()\n",
        "  df_X.replace(cleanup_nums, inplace=True)  # convert continous categorial variables to numerical\n",
        "  df_Y = df[target_field].copy()\n",
        "\n",
        "  x = df_X.to_numpy().astype(np.float32)\n",
        "  y = df_Y.to_numpy().astype(np.float32)\n",
        "\n",
        "  if price_max>0:\n",
        "    idxs = y<price_max\n",
        "    x = x[idxs]\n",
        "    y = y[idxs]\n",
        "\n",
        "  if area_max>0:\n",
        "    idxs = x[:,0]<area_max\n",
        "    x = x[idxs]\n",
        "    y = y[idxs]\n",
        "\n",
        "  return (x, y, df) if return_df else ((x, y, (df_X, df_Y)) if return_df_xy else (x,y))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YqWU0eHts1RM"
      },
      "outputs": [],
      "source": [
        "x, y, df = house_prices_dataset(return_df=True)\n",
        "print(x.shape, y.shape)\n",
        "df.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YDtzVS-1Mxxe"
      },
      "outputs": [],
      "source": [
        "df.describe()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "91nj7znzMEpA"
      },
      "outputs": [],
      "source": [
        "plt.plot(x[:, 0], y, '.')\n",
        "plt.xlabel('area, sq.ft')\n",
        "plt.ylabel('price, $');\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q7CNxkPdNB4L"
      },
      "source": [
        "## 3. Blobs"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "j8wXhleONKgZ"
      },
      "outputs": [],
      "source": [
        "x, y = make_blobs(n_samples=1000, centers=[[0,0], [5,5], [10, 0]])\n",
        "colors = \"ygr\"\n",
        "for i, color in enumerate(colors):\n",
        "    idx = y == i\n",
        "    plt.scatter(x[idx, 0], x[idx, 1], c=color, edgecolor='gray', s=25)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NKcmdcZf0VO8"
      },
      "outputs": [],
      "source": [
        "x, y = make_blobs(n_samples=1000, centers=[[0,0], [5,5], [10, 0]])\n",
        "\n",
        "transformation = [[0.4, 0.2], [-0.4, 1.2]]  # affine transformation matrix\n",
        "x = np.dot(x, transformation)               # applied to point coordinated to make blobs less separable\n",
        "\n",
        "colors = \"ygr\"\n",
        "for i, color in enumerate(colors):\n",
        "    idx = y == i\n",
        "    plt.scatter(x[idx, 0], x[idx, 1], c=color, edgecolor='gray', s=25)\n",
        "plt.show()\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8S1jwU4cXQX4"
      },
      "source": [
        "## 4. MNIST"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e2u82UQ5XQX4"
      },
      "source": [
        "The MNIST database of handwritten digits has a training set of 60,000 examples, and a test set of 10,000 examples. The digits have been size-normalized and centered in a fixed-size image.\n",
        "It is a good database for people who want to try learning techniques and pattern recognition methods on real-world data while spending minimal efforts on preprocessing and formatting (taken from http://yann.lecun.com/exdb/mnist/). Each example is a 28x28 grayscale image and the dataset can be readily downloaded from Tensorflow."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JaNaGGOkXQX5"
      },
      "outputs": [],
      "source": [
        "mnist = tf.keras.datasets.mnist\n",
        "(train_images, train_labels), (test_images, test_labels) = mnist.load_data()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dlUY5gl8XQX7"
      },
      "source": [
        "Let's check few samples:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qtYtGEDdXQX8"
      },
      "outputs": [],
      "source": [
        "n = 3\n",
        "fig, ax = plt.subplots(n, n, figsize=(2*n, 2*n))\n",
        "ax = [ax_xy for ax_y in ax for ax_xy in ax_y]\n",
        "for axi, im_idx in zip(ax, np.random.choice(len(train_images), n**2)):\n",
        "  im = train_images[im_idx]\n",
        "  im_class = train_labels[im_idx]\n",
        "  axi.imshow(im, cmap='gray')\n",
        "  axi.text(1, 4, f'{im_class}', color='r', size=16)\n",
        "  axi.grid(False)\n",
        "plt.tight_layout(pad=0, h_pad=0, w_pad=0)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ITfbaOgfYNsq"
      },
      "source": [
        "## 5. Fashion MNIST"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jgzzOS7YYTru"
      },
      "source": [
        "`Fashion-MNIST` is a dataset of Zalando's article images—consisting of a training set of 60,000 examples and a test set of 10,000 examples. Each example is a 28x28 grayscale image, associated with a label from 10 classes. (from https://github.com/zalandoresearch/fashion-mnist)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RcV2gzmuYljJ"
      },
      "outputs": [],
      "source": [
        "fashion_mnist = tf.keras.datasets.fashion_mnist\n",
        "(train_images, train_labels), (test_images, test_labels) = fashion_mnist.load_data()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SPw6-GoPbT6U"
      },
      "source": [
        "Let's check few samples:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tHFd0sFHY4Li"
      },
      "outputs": [],
      "source": [
        "n = 3\n",
        "fig, ax = plt.subplots(n, n, figsize=(2*n, 2*n))\n",
        "ax = [ax_xy for ax_y in ax for ax_xy in ax_y]\n",
        "for axi, im_idx in zip(ax, np.random.choice(len(train_images), n**2)):\n",
        "  im = train_images[im_idx]\n",
        "  im_class = train_labels[im_idx]\n",
        "  axi.imshow(im, cmap='gray')\n",
        "  axi.text(1, 4, f'{im_class}', color='r', size=16)\n",
        "  axi.grid(False)\n",
        "plt.tight_layout(0,0,0)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "b2LkoWfZEi4g"
      },
      "outputs": [],
      "source": [
        "fmnist_class_names = ['T-shirt/top', 'Trouser', 'Pullover', 'Dress', 'Coat',\n",
        "               'Sandal', 'Shirt', 'Sneaker', 'Bag', 'Ankle boot']"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iHEA0tCLagoV"
      },
      "source": [
        "Each of the training and test examples is assigned to one of the following labels:\n",
        "\n",
        "| Label | Description |\n",
        "| --- | --- |\n",
        "| 0 | T-shirt/top |\n",
        "| 1 | Trouser |\n",
        "| 2 | Pullover |\n",
        "| 3 | Dress |\n",
        "| 4 | Coat |\n",
        "| 5 | Sandal |\n",
        "| 6 | Shirt |\n",
        "| 7 | Sneaker |\n",
        "| 8 | Bag |\n",
        "| 9 | Ankle boot |"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "x2NWxK0BFwyw"
      },
      "source": [
        "## 6. Weather dataset"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LKsTkrMkGB7d"
      },
      "source": [
        "This is a weather time series [dataset](https://www.bgc-jena.mpg.de/wetter/) recorded by the Max Planck Institute for Biogeochemistry\n",
        "It contains weather reacord for 8 years of observation."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "D_bXLpwxGTbZ"
      },
      "outputs": [],
      "source": [
        "def get_weather_df():\n",
        "  # inspired by https://www.tensorflow.org/tutorials/structured_data/time_series\n",
        "\n",
        "  # download and extract dataset\n",
        "  zip_path = tf.keras.utils.get_file(\n",
        "    origin='https://storage.googleapis.com/tensorflow/tf-keras-datasets/jena_climate_2009_2016.csv.zip',\n",
        "    fname='jena_climate_2009_2016.csv.zip',\n",
        "    extract=True)\n",
        "  csv_path, _ = os.path.splitext(zip_path)\n",
        "\n",
        "  # load into pandas df\n",
        "  df = pd.read_csv(csv_path)\n",
        "\n",
        "  # dataset contains records every 10 min, we use hourly records only, thus\n",
        "  # slice [start:stop:step], starting from index 5 take every 6th record\n",
        "  df = df[5::6]\n",
        "\n",
        "  # replace errors in wind velocity to 0\n",
        "  wv = df['wv (m/s)']\n",
        "  bad_wv = wv == -9999.0\n",
        "  wv[bad_wv] = 0.0\n",
        "\n",
        "  max_wv = df['max. wv (m/s)']\n",
        "  bad_max_wv = max_wv == -9999.0\n",
        "  max_wv[bad_max_wv] = 0.0\n",
        "\n",
        "  # obtain timestamps from text time format\n",
        "  date_time = pd.to_datetime(df.pop('Date Time'), format='%d.%m.%Y %H:%M:%S')\n",
        "  timestamp_s = date_time.map(pd.Timestamp.timestamp)\n",
        "  # genarate cyclic features for year and day\n",
        "  day = 24*60*60\n",
        "  year = (365.2425) * day\n",
        "\n",
        "  df['Day sin'] = np.sin(timestamp_s * (2 * np.pi / day))\n",
        "  df['Day cos'] = np.cos(timestamp_s * (2 * np.pi / day))\n",
        "  df['Year sin'] = np.sin(timestamp_s * (2 * np.pi / year))\n",
        "  df['Year cos'] = np.cos(timestamp_s * (2 * np.pi / year))\n",
        "\n",
        "  return df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DsYDOZZKI7_P"
      },
      "outputs": [],
      "source": [
        "weather_df = get_weather_df()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yQHfhrBMJA5w"
      },
      "outputs": [],
      "source": [
        "weather_df.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "wPUvJcY3JjJa"
      },
      "outputs": [],
      "source": [
        "weather_df.describe().T"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "T7qIhgZuMGW4"
      },
      "outputs": [],
      "source": [
        "plt.plot(weather_df['Year cos'])\n",
        "plt.plot(weather_df['Year sin'])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "WxziDnFwLqL0"
      },
      "outputs": [],
      "source": [
        "plt.plot(weather_df['Day cos'][:7*24])\n",
        "plt.plot(weather_df['Day sin'][:7*24])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tDQkIP4cJxuW"
      },
      "outputs": [],
      "source": [
        "def gen_future_T_dataset(X_len=24, Y_offset=48,\n",
        "                         X_features=['p (mbar)', 'T (degC)', 'rh (%)', 'wv (m/s)', 'wd (deg)', 'Day sin', 'Day cos', 'Year sin', 'Year cos'],\n",
        "                         Y_features='T (degC)',\n",
        "                         standardize=True,\n",
        "                         oversample=10\n",
        "                         ):\n",
        "  \"\"\"\n",
        "  Generates pairs of input-label, using sequence of `X_len` samples as input\n",
        "  and value at offset `Y_offset` from start of this sequence as label.\n",
        "  Sample sequnces arte taken at random positions throughout the dataset.\n",
        "  Number of samples is obtained assuming non-overlaping wondows.\n",
        "  Oversampling factor allows to increase this value.\n",
        "\n",
        "  Args:\n",
        "    X_len (int): length of sample sequence\n",
        "    Y_offset (int): offset to the target value from the sequence start\n",
        "    X_features (list or str): features to be used as input\n",
        "    Y_features (list or str): features to be used as labels\n",
        "    standardize (Bool): flag whether to standardize the columns\n",
        "    oversample (int): increases number of samples by this factor wrt baseline\n",
        "                      n = len(df) // (Y_offset+1)\n",
        "  \"\"\"\n",
        "  weather_df = get_weather_df()\n",
        "\n",
        "  if standardize:\n",
        "    mean = weather_df.mean()\n",
        "    std = weather_df.std()\n",
        "    weather_df = (weather_df - mean) / std\n",
        "\n",
        "  df_X = weather_df[X_features]\n",
        "  df_Y = weather_df[Y_features]\n",
        "\n",
        "  n_records = len(df_X)\n",
        "  sample_len = Y_offset+1\n",
        "  n_samples = int((n_records-sample_len)//sample_len*oversample)\n",
        "  offsets = np.random.randint(0, n_records-sample_len, size=n_samples)\n",
        "  offsets.sort()\n",
        "\n",
        "  X = []\n",
        "  Y = []\n",
        "  for o in offsets:\n",
        "    X.append(np.array(df_X[o:o+X_len]))\n",
        "    Y.append(np.array(df_Y[o+Y_offset:o+Y_offset+1]))\n",
        "\n",
        "  X = np.stack(X)\n",
        "  Y = np.concatenate(Y)\n",
        "\n",
        "  return X, Y"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "m-uVf1AOPLIT"
      },
      "outputs": [],
      "source": [
        "x, y = gen_future_T_dataset()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "m2f1pWLfPbJG"
      },
      "outputs": [],
      "source": [
        "x.shape, y.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JJU5MRbdRDEI"
      },
      "outputs": [],
      "source": [
        "for f,fn in enumerate(['p (mbar)', 'T (degC)', 'rh (%)', 'wv (m/s)', 'wd (deg)', 'Day sin', 'Day cos', 'Year sin', 'Year cos']):\n",
        "  plt.figure(figsize=(5*(1+(f==1)), 4))\n",
        "  for s in range(10):\n",
        "    plt.plot(x[s, :, f])\n",
        "    if f==1:\n",
        "      plt.scatter(48, y[s], color=plt.gca().lines[-1].get_color())\n",
        "  plt.title(fn)\n",
        "  plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oAAjJuenj1u0"
      },
      "source": [
        "# 2. Neural Networks Introduction"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5AXSFimKkt91"
      },
      "source": [
        "## 1. Perceptron"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xoFa0c0Xj_6m"
      },
      "source": [
        "(Artificial) Neural network consists of layers of neurons. Artificial neuron, or perceptron, is in fact inspired by a biological neuron.\n",
        "\n",
        "<img src=\"https://github.com/neworldemancer/DSF5/raw/master/figures/Perceptron.png\" alt=\"drawing\" width=\"30%\"/>\n",
        "\n",
        "Such neuron first calculates the linear transformation of the input vector $\\bar x$:\n",
        "$$z = \\bar W \\cdot \\bar x + b = \\sum {W_i x_i} + b$$ where $\\bar W$ is vector of weights and $b$ - bias.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O7cZipYkZlH9"
      },
      "source": [
        "This is effectively linear regression. You can combine those parallely to make a multidimensional prediction, e.g.:\n",
        "\n",
        "<img src=\"https://github.com/neworldemancer/DSF5/raw/master/figures/2Perceptrons.png\" alt=\"drawing\" width=\"30%\"/>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C-mH0li3kzNi"
      },
      "source": [
        "## 2. Nonlinearity"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i74xUp5qkxCy"
      },
      "source": [
        "Combining multiple such objects performing linear transformation sequentially would not bring any additional benefit, as the combined output would still be a linear combination of the inputs.\n",
        "\n",
        "What gives actual power to neurons, is that they additionally perform the nonlinear transformation of the result using activation function $f$ $$y = f(z)$$\n",
        "\n",
        "The most commonly used non-linear transformations are:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iUo59ubLc_fi"
      },
      "outputs": [],
      "source": [
        "def ReLU(z):\n",
        "  return np.clip(z, a_min=0, a_max=np.max(z))\n",
        "def SELU(z, a=1):\n",
        "  p = np.clip(z, a_min=0, a_max=np.max(z))\n",
        "  n = np.clip(z, a_min=np.min(z), a_max=0)\n",
        "  return p + (np.exp(n)-1) * a\n",
        "def LReLU(z, a=0.1):\n",
        "  return np.clip(z, a_min=0, a_max=np.max(z)) + np.clip(z, a_min=np.min(z), a_max=0) * a\n",
        "def sigmoid(z):\n",
        "  return 1/(1 + np.exp(-z))\n",
        "def step(z):\n",
        "  return np.heaviside(z, 0)\n",
        "fig, ax = plt.subplots(1, 6, figsize=(18, 3))\n",
        "z = np.linspace(-10, 10, 100)\n",
        "ax[0].plot(z, ReLU(z))\n",
        "ax[0].set_title('Rectified Linear Unit (LU)')\n",
        "ax[1].plot(z, LReLU(z))\n",
        "ax[1].set_title('Leaky Rectified LU')\n",
        "ax[2].plot(z, SELU(z))\n",
        "ax[2].set_title('Scaled Exponential LU')\n",
        "ax[3].plot(z, sigmoid(z))\n",
        "ax[3].set_title(r'$\\sigma$(z)=$\\frac{1}{1+e^z}$')\n",
        "ax[4].plot(z, np.tanh(z))\n",
        "ax[4].set_title('Hyperbolic tangent');\n",
        "ax[5].plot(z, step(z))\n",
        "ax[5].text(-6, 0.5, 'NOT USED', size=19, c='r')\n",
        "ax[5].set_title('Step function');\n",
        "for axi in ax:\n",
        "  axi.set_xlabel('z')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QsEbaIujf0DV"
      },
      "source": [
        "And the reason we don't use a simple step function, is that it's everywhere either not differentiable or its derivative is zero.\n",
        "\n",
        "The last nonlinearity to mention here is *softmax*:\n",
        "$$y_i = SoftMax(\\bar z)_i = \\frac{ e^{z_i}}{\\sum_j e^{z_j}}$$\n",
        "\n",
        "While each $z_i$ can have any value, the corresponding $y_i\\in[0,1]$, and $\\sum_i y_i=1$, just like probabilities!\n",
        "\n",
        "While these $y_i$ are only pseudo-probabilities, this nonlinearity allows one to model probabilities, e.g. of a data-point belonging to a certain class.\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8QZtmo6Vk2rp"
      },
      "source": [
        "## 3. Fully connected net"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PM5z2czOc_4A"
      },
      "source": [
        "In a fully connected neural network each layer is a set of N neurons, performing different transformations of all the same layer's inputs $\\bar x = [x_i]$ producing output vector $\\bar y = [y_j]_{i=1..N}$: $$y_j = f(\\bar W_j \\cdot \\bar x + b_j)$$\n",
        "\n",
        "Since output of each layer forms input of next layer, one can write for layer $l$ (upper index denotes layer): $$x^l_j = f(\\bar W^l_j \\cdot \\bar x^{l-1} + b^l_j)$$ where $\\bar x^0$ is network's input vector.\n",
        "\n",
        "<img src=\"https://github.com/neworldemancer/DSF5/raw/master/figures/MLP.png\" alt=\"drawing\" width=\"50%\"/>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u62hCoFbklaW"
      },
      "source": [
        "## 4. Loss function"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-58_mV6ElY5C"
      },
      "source": [
        "The last part of the puzzle is the measure of network performance, which is used to optimize the network's parameters $W^l_j$ and $b^l_j$.\n",
        "Denoting the network's output for an input $x_i$ as $\\hat y_i=\\hat y_i(x_i)$ and given the label $y_i$:\n",
        "\n",
        "1. In case of regression loss shows \"distance\" from target values:\n",
        "* L2 (MSE): $L = \\sum_i (y_i-\\hat y_i)^2$\n",
        "* L1 (MAE): $L = \\sum_i |y_i-\\hat y_i|$\n",
        "\n",
        "1. In case of classification we can use cross-entropy, which shows \"distance\" from target distribution:\n",
        "$$L = - \\sum_i \\sum_c y_{i,c} \\log(\\hat y_{i,c})$$\n",
        "Here $\\hat y_{i,c}$ - pseudo-probability of $x_i$ belonging to class $c$ and $y_{i,c}$ uses 1-hot encoding:\n",
        "\n",
        "$$y_{i,c}=\n",
        "\\begin{cases}\n",
        "    1,& \\text{if } x_i \\text{ belongs to class } c\\\\\n",
        "    0,              & \\text{otherwise}\n",
        "\\end{cases}$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "--9TcJQ8Rahm"
      },
      "source": [
        "## 5. Training"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "urF63xIfRah0"
      },
      "source": [
        "Training of neural networks is performed iteratively. The weights  $W^l_j$ and $b^l_j$ are updated on each iteration of training according to the value of the derivative of the loss function with respect to corresponding parameter:\n",
        "$$W^l_j \\rightarrow W^l_j - \\lambda \\frac{\\partial L}{\\partial W^l_j }$$\n",
        "$$b^l_j \\rightarrow b^l_j - \\lambda \\frac{\\partial L}{\\partial b^l_j },$$\n",
        "\n",
        "This is Gradient Descent optimization with learning rate $\\lambda$. The partial derivatives are calculated by the chain law, and this approach is known as [backpropagation](https://en.wikipedia.org/wiki/Backpropagation).\n",
        "\n",
        "In practice often for each iteration the loss $L$ is evaluated not on all samples, but on a sub-sample, so-called *minibatch* (sometimes - just batch). In most cases the sample order and selection for each minibatch is performed at random rendering this approach to be stochastic (thus it's called [Stochastic Gradient Descent](https://en.wikipedia.org/wiki/Stochastic_gradient_descent)). One iteration through all training data in minibatches is called *epoch*."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yrwaonuB4F8m"
      },
      "source": [
        "# 3. Regression with neural network"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Yj48cHUiVttY"
      },
      "source": [
        "## 2. 1D input"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_-SUw6f7Vwo7"
      },
      "outputs": [],
      "source": [
        "x = np.linspace(0, 10, num=1000)\n",
        "y = ((x>3)*(x<5)).astype(np.float32)\n",
        "\n",
        "#y = x**3/500 + np.sin(x) #+ np.sin((10-x)**2)\n",
        "\n",
        "x /= 10\n",
        "x -= 0.5\n",
        "plt.plot(x, y);"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "M3nmmED5WRo0"
      },
      "outputs": [],
      "source": [
        "model = tf.keras.models.Sequential([\n",
        "  tf.keras.layers.Flatten(input_shape=(1,)),\n",
        "  tf.keras.layers.Dense(100, activation='relu'),\n",
        "  tf.keras.layers.Dense(1, activation=None),\n",
        "])\n",
        "model.compile(optimizer=tf.keras.optimizers.Adam(learning_rate=0.01),\n",
        "              loss='mae',\n",
        "              metrics=['mse'])\n",
        "model.summary()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Z3RkZ057yEo8"
      },
      "outputs": [],
      "source": [
        "class saving_pred(tf.keras.callbacks.Callback):\n",
        "  def __init__(self):\n",
        "    self.y_preds = []\n",
        "\n",
        "  def on_epoch_end(self, epoch, logs=None):\n",
        "    if epoch % 10 == 0:\n",
        "        y_p = model.predict(x).flatten()\n",
        "        #print('\\n', epoch, y_p.shape, '\\n\\n')\n",
        "        self.y_preds.append(y_p)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3Gw8BJYxWWql"
      },
      "outputs": [],
      "source": [
        "saving_callback = saving_pred()\n",
        "hist = model.fit(x, y,\n",
        "                 epochs=300, batch_size=100, callbacks=[saving_callback])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yLd8yWdZoJD7"
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(1, 2, figsize=(10,5))\n",
        "axs[0].plot(hist.epoch, hist.history['loss'])\n",
        "axs[0].set_title('loss')\n",
        "axs[1].plot(hist.epoch, hist.history['mse'])\n",
        "axs[1].set_title('mse')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_UhGnieHWcP8"
      },
      "outputs": [],
      "source": [
        "y_p = model.predict(x)\n",
        "im = plt.plot(x, y_p)\n",
        "im = plt.plot(x, y, '--')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "t4Lh2YwDybnl"
      },
      "outputs": [],
      "source": [
        "\n",
        "%%capture\n",
        "y_preds = saving_callback.y_preds\n",
        "\n",
        "fig = plt.figure()\n",
        "line = plt.plot(x, y_preds[0])\n",
        "line2 = plt.plot(x, y)\n",
        "\n",
        "def animate(i):\n",
        "    y_p = y_preds[i]\n",
        "    line[0].set_data((x, y_p))\n",
        "    return line\n",
        "\n",
        "ani = animation.FuncAnimation(fig, animate, frames=len(y_preds))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0DTtFVuDr3b8"
      },
      "outputs": [],
      "source": [
        "ani"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2PFZcbkbVprF"
      },
      "source": [
        "## 2. 2D input"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zUCua55tE85o"
      },
      "source": [
        "Here we will build a neural network to perform regression in 2D: we will fit $(x,y)$ pixel coordinates of an image to the pixel brightness at that location."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "a6rG4pXiwT1o"
      },
      "outputs": [],
      "source": [
        "url = 'https://github.com/neworldemancer/DSF5/raw/master/figures/unibe.jpg'\n",
        "image_big = imread(url)\n",
        "image_big = image_big[...,0:3]/255\n",
        "plt.imshow(image_big)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3MzVfGJocmDo"
      },
      "outputs": [],
      "source": [
        "def subsample(im, factor):\n",
        "  h, w = im.shape[:2]\n",
        "  h = int(h/factor)\n",
        "  w = int(w/factor)\n",
        "\n",
        "  img = Image.fromarray((im*255).astype('uint8'))\n",
        "  newimg = img.resize((w, h), Image.LANCZOS)\n",
        "  return np.asarray(newimg)/255"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9iGxdgfwEeiD"
      },
      "outputs": [],
      "source": [
        "image = subsample(image_big, 15)\n",
        "image = image.mean(axis=2, keepdims=True)\n",
        "plt.imshow(image[...,0], cmap='gray')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9hGEsYDX4hwu"
      },
      "outputs": [],
      "source": [
        "h, w, c = image.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0RiEARRw4wAP"
      },
      "outputs": [],
      "source": [
        "X0 = np.meshgrid(np.linspace(0, 1, w), np.linspace(0, 1, h))\n",
        "X = np.stack(X0, axis=-1).reshape((-1, 2))\n",
        "\n",
        "Y = image.reshape((-1, c))\n",
        "X.shape, Y.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "sBVWYpq75GSa"
      },
      "outputs": [],
      "source": [
        "model = tf.keras.models.Sequential([\n",
        "  tf.keras.layers.Flatten(input_shape=(2,)),\n",
        "  tf.keras.layers.Dense(c, activation='sigmoid'),\n",
        "])\n",
        "model.compile(optimizer='adam',\n",
        "              loss='mae',\n",
        "              metrics=['mse'])\n",
        "model.summary()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tR7OlLN2XEo1"
      },
      "outputs": [],
      "source": [
        "hist = model.fit(X, Y, epochs=500, batch_size=2048)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7xUI4w3zXdoy"
      },
      "outputs": [],
      "source": [
        "Y_p = model.predict(X)\n",
        "Y_p = Y_p.reshape((h,w,c))\n",
        "im = plt.imshow(Y_p[...,0], cmap='gray')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "i9yeAZ9zaQH9"
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(1, 2, figsize=(10,5))\n",
        "axs[0].plot(hist.epoch, hist.history['loss'])\n",
        "axs[0].set_title('loss')\n",
        "axs[1].plot(hist.epoch, hist.history['mse'])\n",
        "axs[1].set_title('mse')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sxw90MRNDJu9"
      },
      "source": [
        "What is actually happening here? We fit an image with an $\\mathbb{R}^2 \\rightarrow \\mathbb{R}$ funcion."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X27VLMYS8R8D"
      },
      "outputs": [],
      "source": [
        "\n",
        "x = X[:,0].reshape((h,w))\n",
        "y = X[:,1].reshape((h,w))\n",
        "z = Y.reshape((h,w))\n",
        "zp = Y_p.reshape((h,w))\n",
        "\n",
        "ds = 2\n",
        "x = subsample(x, ds)\n",
        "y =-subsample(y, ds)\n",
        "z = subsample(z, ds)\n",
        "zp = subsample(zp, ds)\n",
        "\n",
        "#surf = ax.plot_surface(x,y,z, cmap='coolwarm', linewidth=0, antialiased=False, alpha=0.3)\n",
        "#surf = ax.plot_surface(x,y,zp, cmap='coolwarm', linewidth=0, antialiased=False, alpha=0.5)\n",
        "#plt.show()\n",
        "\n",
        "fig = go.Figure(data=[go.Surface(z=zp, x=x, y=y)])\n",
        "fig.update_layout(autosize=False,\n",
        "                  width=500, height=500,\n",
        "                  margin=dict(l=65, r=50, b=65, t=90))\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ONLyx-kYXwQe"
      },
      "source": [
        "Let's try the same with an RGB image:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "WR3Te8_NX43W"
      },
      "outputs": [],
      "source": [
        "image = subsample(image_big, 15)\n",
        "plt.imshow(image)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gubPoLVcbzdB"
      },
      "outputs": [],
      "source": [
        "h, w, c = image.shape\n",
        "X = np.meshgrid(np.linspace(0, 1, w), np.linspace(0, 1, h))\n",
        "X = np.stack(X, axis=-1).reshape((-1, 2))\n",
        "\n",
        "Y = image.reshape((-1, c))\n",
        "X.shape, Y.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Kg1iKZ8LXl1R"
      },
      "outputs": [],
      "source": [
        "model = tf.keras.models.Sequential([\n",
        "  tf.keras.layers.Flatten(input_shape=(2,)),\n",
        "  tf.keras.layers.Dense(c, activation='sigmoid'),\n",
        "])\n",
        "model.compile(optimizer='adam',\n",
        "              loss='mae',\n",
        "              metrics=['mse'])\n",
        "model.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YkLrIjf0YAUF"
      },
      "source": [
        "But now we will save images during the course of training, at first after 2, then 4, 8, 16, etc epochs.\n",
        "(**Remember**: call to `model.fit` does NOT reinitialize trainable variables. Every time it continues from the previous state):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4aelmQDvzunH"
      },
      "outputs": [],
      "source": [
        "class saving_pred(tf.keras.callbacks.Callback):\n",
        "  def __init__(self, n_ep_pow):\n",
        "    self.ims = []\n",
        "    self.save_epochs = [0]+[2**i for i in range(n_ep_pow)]\n",
        "\n",
        "  def on_epoch_end(self, epoch, logs=None):\n",
        "    if epoch in self.save_epochs:\n",
        "        Y_p = model.predict(X)\n",
        "        Y_p = Y_p.reshape((h, w, c))\n",
        "        self.ims.append(Y_p)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0ezFwya85m4M"
      },
      "outputs": [],
      "source": [
        "n_ep_pow = 11\n",
        "n_ep = 2**n_ep_pow+1\n",
        "\n",
        "ims_callback = saving_pred(n_ep_pow)\n",
        "hist = model.fit(X, Y, epochs=n_ep, batch_size=1*2048, verbose=0, callbacks=[ims_callback])\n",
        "\n",
        "plt.plot(hist.epoch, hist.history['loss'])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xFdyqtEyx3-f"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "\n",
        "ims = ims_callback.ims\n",
        "fig = plt.figure()\n",
        "im = plt.imshow(ims[0])\n",
        "\n",
        "def animate(i):\n",
        "    img = ims[i]\n",
        "    im.set_data(img)\n",
        "    return im\n",
        "\n",
        "ani = animation.FuncAnimation(fig, animate, frames=len(ims))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TxFYMEy5WmF4"
      },
      "outputs": [],
      "source": [
        "ani"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1wc9tXb9a8_q"
      },
      "source": [
        "While the colors properly represent the target image, our model still poses very limited capacity, allowing it to effectively represent only 3 boundaries.\n",
        "\n",
        "Let's upscale our model:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mHmrapmHbWxN"
      },
      "outputs": [],
      "source": [
        "model = tf.keras.models.Sequential([\n",
        "  tf.keras.layers.Flatten(input_shape=(2,)),\n",
        "  tf.keras.layers.Dense(128, activation='relu'),\n",
        "  tf.keras.layers.Dense(8, activation='relu'),\n",
        "  tf.keras.layers.Dense(c, activation='sigmoid'),\n",
        "])\n",
        "model.compile(optimizer=tf.keras.optimizers.Adam(learning_rate=0.005), # optimizer='adam',\n",
        "              loss='mae',\n",
        "              metrics=['mse'])\n",
        "model.summary()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-uqJWwiu2B9r"
      },
      "outputs": [],
      "source": [
        "X.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GmAqFs7MfThU"
      },
      "outputs": [],
      "source": [
        "n_ep_pow = 13\n",
        "n_ep = 2**n_ep_pow+1\n",
        "\n",
        "ims_callback = saving_pred(n_ep_pow)\n",
        "hist = model.fit(X, Y, epochs=n_ep, batch_size=1*2048, verbose=0, callbacks=[ims_callback])\n",
        "\n",
        "plt.plot(hist.epoch, hist.history['loss'])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7GruTrQNeC4A"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "\n",
        "ims = ims_callback.ims\n",
        "fig = plt.figure()\n",
        "im = plt.imshow(ims[0])\n",
        "\n",
        "def animate(i):\n",
        "    img = ims[i]\n",
        "    im.set_data(img)\n",
        "    return im\n",
        "\n",
        "ani = animation.FuncAnimation(fig, animate, frames=len(ims))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mCAOLpsMeNRA"
      },
      "outputs": [],
      "source": [
        "ani"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rEVkc53u3Vz_"
      },
      "source": [
        "We can also now render the fitted image at any resolution, by evaluating values in between of the original pixels:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yV0q97Ip3dik"
      },
      "outputs": [],
      "source": [
        "scale = 30\n",
        "w_rend = w * scale\n",
        "h_rend = h * scale\n",
        "X_rend = np.meshgrid(np.linspace(0, 1, w_rend), np.linspace(0, 1, h_rend))\n",
        "X_rend = np.stack(X_rend, axis=-1).reshape((-1, 2))\n",
        "\n",
        "Y_rend = model.predict(X_rend)\n",
        "Y_rend = Y_rend.reshape((h_rend, w_rend, c))\n",
        "\n",
        "plt.imshow(Y_rend)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vt1q-Nez5Hak"
      },
      "outputs": [],
      "source": [
        "plt.imshow(image_big)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pCFB8rQ15Tf9"
      },
      "outputs": [],
      "source": [
        "plt.imshow(image)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aOk-BLOFeV61"
      },
      "source": [
        "## EXERCISE 1. Regression on an image with Neural Network."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "09yEvb4aebd5"
      },
      "source": [
        "Load some image, downscale to a similar resolution, and train a deeper model, for example 5 layers, more parameters in widest layers."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "K415Bvdguld7"
      },
      "outputs": [],
      "source": [
        "# 1. Load your image\n",
        "\n",
        "# 2. build a deeper model\n",
        "\n",
        "# 3. inspect the evolution"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FCQ5wxQCgVXD"
      },
      "source": [
        "# 4. Classification of the F-MNIST dataset with neural network"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dSRc4lzeYgxX"
      },
      "source": [
        "## 0. Load data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xM_660aCYmnS"
      },
      "source": [
        "We will create a model for classification of the F-MNIST dataset that we go acquainted with in previous sessions. We will normalize the inputs to have values $\\in[0,1]$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "maLerZKejJxG"
      },
      "outputs": [],
      "source": [
        "fashion_mnist = tf.keras.datasets.fashion_mnist\n",
        "(x_train, y_train), (x_test, y_test) = fashion_mnist.load_data()\n",
        "x_train = x_train/255\n",
        "x_test = x_test/255\n",
        "\n",
        "class_names = ['T-shirt/top', 'Trouser', 'Pullover', 'Dress', 'Coat',\n",
        "               'Sandal', 'Shirt', 'Sneaker', 'Bag', 'Ankle boot']"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "agJ3gkg0tkTJ"
      },
      "source": [
        "## 1. Building a neural network"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "D0dtcC9MZG3N"
      },
      "outputs": [],
      "source": [
        "print(x_train[0].shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cjVnBkWjZHVf"
      },
      "source": [
        "The size of each image sample $-\\; 28\\times28\\text{ pixels}\\;-\\;$ defines the input size for our neural network. Network's output - probabilities of belonging to each of the 10 classes."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lPs6sljtjDUE"
      },
      "source": [
        "The following creates a 'model'. It is an object containing the neural network model itself - a simple 3-layer fully connected neural network, optimization parameters, as well as the interface for model training."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "code",
        "id": "LC1SpGLbtkTL"
      },
      "outputs": [],
      "source": [
        "model = tf.keras.models.Sequential([\n",
        "  tf.keras.layers.Flatten(input_shape=(28, 28)),\n",
        "  tf.keras.layers.Dense(10, activation='softmax')\n",
        "])\n",
        "\n",
        "model.compile(optimizer='adam',\n",
        "              loss='sparse_categorical_crossentropy',\n",
        "              metrics=['accuracy'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RFwr5MLMjxI1"
      },
      "source": [
        "Model summary provides information about the model's layers and trainable parameters"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bCttp5zeb5l2"
      },
      "outputs": [],
      "source": [
        "model.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "P18eyAQHqZGG"
      },
      "source": [
        "## 2. Model training"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VNIdb5Gtlr32"
      },
      "source": [
        "The `fit` function is the interface for model training.\n",
        "Here one can specify training and validation datasets, minibatch size, and the number of training epochs.\n",
        "\n",
        "Here during training we also save the trained models checkpoints after each epoch of training."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0-OxT0aNVx-y"
      },
      "outputs": [],
      "source": [
        "save_path = 'save/mnist_{epoch}.weights.h5'\n",
        "save_callback = tf.keras.callbacks.ModelCheckpoint(filepath=save_path, save_weights_only=True)\n",
        "\n",
        "hist = model.fit(x=x_train, y=y_train,\n",
        "                 epochs=50, batch_size=128,\n",
        "                 validation_data=(x_test, y_test),\n",
        "                 callbacks=[save_callback])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8l9Gz1e4V-7Q"
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(1, 2, figsize=(10,5))\n",
        "axs[0].plot(hist.epoch, hist.history['loss'])\n",
        "axs[0].plot(hist.epoch, hist.history['val_loss'])\n",
        "axs[0].legend(('training loss', 'validation loss'), loc='lower right')\n",
        "axs[1].plot(hist.epoch, hist.history['accuracy'])\n",
        "axs[1].plot(hist.epoch, hist.history['val_accuracy'])\n",
        "\n",
        "axs[1].legend(('training accuracy', 'validation accuracy'), loc='lower right')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fUnlYrfaBmQ8"
      },
      "source": [
        "Current model performance can be evaluated, e.g. on the test dataset:\n",
        "\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_cq_gqG4V9il"
      },
      "outputs": [],
      "source": [
        "model.evaluate(x_test,  y_test, verbose=2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4-5qgb0rDyj4"
      },
      "source": [
        "We can test trained model on an image:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "BPU8mOg2DVSO"
      },
      "outputs": [],
      "source": [
        "im_id = 0\n",
        "y_pred = model(x_test)\n",
        "\n",
        "y_pred_most_probable = np.argmax(y_pred[im_id])\n",
        "print('true lablel: ', y_test[im_id],\n",
        "      '; predicted: ',  y_pred_most_probable,\n",
        "      f'({class_names[y_pred_most_probable]})')\n",
        "plt.imshow(x_test[im_id], cmap='gray');"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3uPb0WPTk6oq"
      },
      "source": [
        "As well as inspect on which samples does the model fail:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OKx-umE7R6AL"
      },
      "outputs": [],
      "source": [
        "y_pred_most_probable_all = np.argmax(y_pred, axis=1)\n",
        "wrong_pred_map = y_pred_most_probable_all!=y_test\n",
        "wrong_pred_idx = np.arange(len(wrong_pred_map))[wrong_pred_map]\n",
        "\n",
        "im_id = wrong_pred_idx[0]\n",
        "\n",
        "y_pred_most_probable = y_pred_most_probable_all[im_id]\n",
        "print('true lablel: ', y_test[im_id],\n",
        "      f'({class_names[y_test[im_id]]})',\n",
        "      '; predicted: ',  y_pred_most_probable,\n",
        "      f'({class_names[y_pred_most_probable]})')\n",
        "plt.imshow(x_test[im_id], cmap='gray');\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HkrIXxYmqiyR"
      },
      "source": [
        "## 3. Loading trained model"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gt_BelVEdH1y"
      },
      "outputs": [],
      "source": [
        "model.load_weights('save/mnist_1.weights.h5')\n",
        "model.evaluate(x_test,  y_test, verbose=2)\n",
        "\n",
        "model.load_weights('save/mnist_12.weights.h5')\n",
        "model.evaluate(x_test,  y_test, verbose=2)\n",
        "\n",
        "model.load_weights('save/mnist_18.weights.h5')\n",
        "model.evaluate(x_test,  y_test, verbose=2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MxIIoNRYqqd_"
      },
      "source": [
        "## 4. Inspecting trained variables"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5waklWBUBwuO"
      },
      "source": [
        "We can obtain the trained variables from model layers:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "i-nnCph8rU01"
      },
      "outputs": [],
      "source": [
        "l = model.get_layer(index=1)\n",
        "w, b = l.weights\n",
        "\n",
        "w = w.numpy()\n",
        "b = b.numpy()\n",
        "print(w.shape, b.shape)\n",
        "w = w.reshape((28,28,-1)).transpose((2, 0, 1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Hfmd5YeUCCKO"
      },
      "source": [
        "Let's visualize first 10:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "q8UOgBWJfzMg"
      },
      "outputs": [],
      "source": [
        "n = 10\n",
        "fig, axs = plt.subplots(1, n, figsize=(4.1*n,4))\n",
        "for i, wi in enumerate(w[:n]):\n",
        "  axs[i].imshow(wi, cmap='gray')\n",
        "  axs[i].set_title(class_names[i])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i02IbJ2gtkTT"
      },
      "source": [
        "## 6. Inspecting gradients"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-uuktPZ9CX8W"
      },
      "source": [
        "We can also evaluate the gradients of each output with respect to an input:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "L371n9COtkTU"
      },
      "outputs": [],
      "source": [
        "idx = 112\n",
        "inp_v = x_train[idx:idx+1]  # use some image to compute gradients with respect to\n",
        "\n",
        "inp = tf.constant(inp_v)  # create tf constant tensor\n",
        "with tf.GradientTape() as tape:  # gradient tape for gradint evaluation\n",
        "  tape.watch(inp)  # take inp as variable\n",
        "  preds = model(inp) # evaluate model output\n",
        "\n",
        "grads = tape.jacobian(preds, inp)  # evaluate d preds[i] / d inp[j]\n",
        "print(grads.shape, '<- (Batch_preds, preds[i], Batch_inp, inp[y], inp[x])')\n",
        "grads = grads.numpy()[0,:,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "83mY_1BKiIIB"
      },
      "outputs": [],
      "source": [
        "print('prediction:', np.argmax(preds[0]))\n",
        "fig, axs = plt.subplots(1, 11, figsize=(4.1*11,4))\n",
        "axs[0].imshow(inp_v[0])\n",
        "axs[0].set_title('raw')\n",
        "vmin,vmax = grads.min(), grads.max()\n",
        "for i, g in enumerate(grads):\n",
        "  axs[i+1].imshow(g, cmap='gray', vmin=vmin, vmax=vmax)\n",
        "  axs[i+1].set_title(r'$\\frac{\\partial\\;P(digit\\,%d)}{\\partial\\;input}$' % i, fontdict={'size':16})\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TDS71kLLmRE2"
      },
      "source": [
        "## EXERCISE 2: Train deeper network for F-MNIST classification"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8eyO5lzOmWAG"
      },
      "source": [
        "Make a deeper model, with wider layers. Remember to use the `'softmax'` activation in the last layer, as required for the classification task to encode pseudoprobabilities. In the other layers you could use `'relu'`.\n",
        "\n",
        "Try to achieve 90% accuracy.\n",
        "Does your model overfit?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8CUw7uHUxSoQ"
      },
      "outputs": [],
      "source": [
        "# 1. create model\n",
        "# 2. train the model\n",
        "# 3. plot the loss and accuracy evolution during training\n",
        "# 4. evaluate model in best point (before overfitting)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CWm4FgSu-tK0"
      },
      "source": [
        "# 5. Extras and Q&A"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OReFhRTHxCdl"
      },
      "source": [
        "<img src=\"https://github.com/neworldemancer/DSF5/raw/master/figures/cheatsheet.png\" width=\"100%\"/>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_65g5ZdTxCdm"
      },
      "source": [
        "<img src=\"https://github.com/neworldemancer/DSF5/raw/master/figures/clusters.png\" width=\"100%\"/>"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "LzEnLAMt3De2",
        "pclZR6uFklf_",
        "8UQgU5I-lEll",
        "FJ5rjq7fIe8Q",
        "q7CNxkPdNB4L",
        "8S1jwU4cXQX4",
        "ITfbaOgfYNsq",
        "x2NWxK0BFwyw",
        "oAAjJuenj1u0",
        "5AXSFimKkt91",
        "C-mH0li3kzNi",
        "8QZtmo6Vk2rp",
        "u62hCoFbklaW",
        "--9TcJQ8Rahm",
        "yrwaonuB4F8m",
        "Yj48cHUiVttY",
        "2PFZcbkbVprF",
        "FCQ5wxQCgVXD",
        "dSRc4lzeYgxX",
        "agJ3gkg0tkTJ",
        "CWm4FgSu-tK0"
      ],
      "provenance": []
    },
    "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.7.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}