tutorial_training.ipynb

{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "# Tutorial: Training a model\n",
    "In this tutorial, we will train an event-based state-space model on a reduced version of the [Spiking Heidelberg Digits](https://zenkelab.org/resources/spiking-heidelberg-datasets-shd/) dataset.\n",
    "For training on larger datasets or multiple GPUs, we recommend using the training script `run_training.py` instead.\n",
    "\n",
    "## Setup\n",
    "\n",
    "Install and load the important modules and configuration. To install required packages, please do \n",
    "```\n",
    "pip3 install requirements.txt\n",
    "```\n",
    "\n",
    "Directories for loading datasets, model checkpoints and saving results are defined in the configuration file `system/local.yaml`.\n",
    "Please set your directories accordingly."
   ],
   "id": "4d02d51dcadfcfb"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Data loading\n",
    "The SHD dataset contains 20 classes, digits from 0 to 9 in both German and English. \n",
    "We will use a reduced version of the dataset containing only two digits to train the model to non-trivial performance in reasonable time even on CPUs.\n",
    "\n",
    "[Download the training and test dataset](https://zenkelab.org/datasets/) and unpack the archives to `./data/`."
   ],
   "id": "df7ee68ff3e429ed"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "from torch.utils.data import Dataset, DataLoader, random_split\n",
    "import h5py\n",
    "import numpy as np\n",
    "\n",
    "class SpikingHeidelbergDigits(Dataset):\n",
    "    def __init__(self, path_to_file):\n",
    "        self.num_classes = 2\n",
    "        self.num_channels = 700\n",
    "        self.path_to_file = path_to_file\n",
    "        \n",
    "        # load the dataset\n",
    "        with h5py.File(path_to_file, 'r') as f:\n",
    "            self.channels = f['spikes']['units'][:]\n",
    "            self.timesteps = f['spikes']['times'][:]\n",
    "            self.labels = f['labels'][:]\n",
    "        \n",
    "        # filter the dataset to contain only two classes\n",
    "        mask = (self.labels == 0) | (self.labels == 1)\n",
    "        self.channels = self.channels[mask]\n",
    "        self.timesteps = self.timesteps[mask]\n",
    "        self.labels = self.labels[mask]\n",
    "        \n",
    "    def __len__(self):\n",
    "        return len(self.labels)\n",
    "    \n",
    "    def __getitem__(self, idx):\n",
    "        # create tonic-like structured arrays\n",
    "        dtype = np.dtype([(\"t\", int), (\"x\", int), (\"p\", int)])\n",
    "        struct_arr = np.empty_like(self.channels[idx], dtype=dtype)\n",
    "        \n",
    "        # yield timesteps in milliseconds\n",
    "        timesteps = self.timesteps[idx] * 1e6\n",
    "        \n",
    "        struct_arr['t'] = timesteps\n",
    "        struct_arr['x'] = self.channels[idx]\n",
    "        struct_arr['p'] = 1\n",
    "        \n",
    "        # one-hot encoding of labels (required for CutMix augmentation)\n",
    "        label = np.eye(self.num_classes)[self.labels[idx]].astype(np.int32)\n",
    "            \n",
    "        return struct_arr, label"
   ],
   "id": "f9883d23c86e5bcd",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Load the training and test dataset\n",
    "train_dataset = SpikingHeidelbergDigits('data/shd_train.h5')\n",
    "test_dataset = SpikingHeidelbergDigits('data/shd_test.h5')"
   ],
   "id": "3be0429979f96a3f",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Check the length of the datasets to check if the data loading was successful.",
   "id": "cf72529578541b9"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "print(f\"Number of training samples: {len(train_dataset)}\")\n",
    "print(f\"Number of test samples: {len(test_dataset)}\")"
   ],
   "id": "ec059aefa5d3408",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Now, create a validation set by randomly splitting the training dataset, and create data loaders for training, validation, and test datasets.",
   "id": "ee746256534411df"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# Split the training dataset into training and validation\n",
    "train_dataset, val_dataset = random_split(train_dataset, [int(0.8*len(train_dataset)), len(train_dataset) - int(0.8*len(train_dataset))])\n",
    "\n",
    "# Create data loaders\n",
    "from event_ssm.dataloading import event_stream_collate_fn\n",
    "from functools import partial\n",
    "\n",
    "collate_fn = partial(event_stream_collate_fn, resolution=(700,), pad_unit=8192)\n",
    "train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True, drop_last=True, collate_fn=collate_fn)\n",
    "val_loader = DataLoader(val_dataset, batch_size=32, shuffle=False, drop_last=True, collate_fn=collate_fn)\n",
    "test_loader = DataLoader(test_dataset, batch_size=32, shuffle=False, collate_fn=collate_fn)"
   ],
   "id": "1ab24a1d63c4c194",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Model definition\n",
    "We use the [hydra](https://hydra.cc/docs/intro/) package for efficient configuration management. Define the model configuration in a config file in the `configs` directory."
   ],
   "id": "acc88e3270fda10b"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "from hydra import compose, initialize\n",
    "from omegaconf import OmegaConf, open_dict\n",
    "\n",
    "with initialize(version_base=None, config_path=\"configs\", job_name=\"training tutorial\"):\n",
    "    cfg = compose(config_name=\"base\", overrides=[\"task=tutorial\"])\n",
    "\n",
    "with open_dict(cfg):    \n",
    "    # optax updates the schedule every iteration and not every epoch\n",
    "    cfg.optimizer.total_steps = cfg.training.num_epochs * len(train_loader) // cfg.optimizer.accumulation_steps\n",
    "    cfg.optimizer.warmup_steps = cfg.optimizer.warmup_epochs * len(train_loader) // cfg.optimizer.accumulation_steps\n",
    "    \n",
    "    # scale learning rate by batch size\n",
    "    cfg.optimizer.ssm_lr = cfg.optimizer.ssm_base_lr * cfg.training.per_device_batch_size * cfg.optimizer.accumulation_steps\n",
    "\n",
    "print(OmegaConf.to_yaml(cfg))"
   ],
   "id": "810edc0798ad7622",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Now, create the model using the configuration defined above.",
   "id": "2ca62d33ebabfdb2"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "from event_ssm.ssm import init_S5SSM\n",
    "from event_ssm.seq_model import BatchClassificationModel\n",
    "\n",
    "ssm_init_fn = init_S5SSM(**cfg.model.ssm_init)\n",
    "model = BatchClassificationModel(\n",
    "    ssm=ssm_init_fn,\n",
    "    num_classes=test_dataset.num_classes,\n",
    "    num_embeddings=test_dataset.num_channels,\n",
    "    **cfg.model.ssm,\n",
    ")"
   ],
   "id": "83e2062ea8b4fe02",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "\n",
    "Initialize the training state by feeding a dummy input"
   ],
   "id": "c855737447f70896"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import jax\n",
    "from event_ssm.train_utils import init_model_state\n",
    "\n",
    "# pick the first batch from the training loader\n",
    "batch = next(iter(train_loader))\n",
    "inputs, targets, timesteps, lengths = batch\n",
    "\n",
    "# initialize the training state\n",
    "key = jax.random.PRNGKey(cfg.seed)\n",
    "state = init_model_state(key, model, inputs, timesteps, lengths, cfg.optimizer)"
   ],
   "id": "d4def1c65952a8ba",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Inspect the model\n",
    "The model parameters are accessible as part of the training state. \n",
    "We will look into the spectrum of the recurrent operator here.\n",
    "The model was initialized with a single stage of blocks."
   ],
   "id": "424bce6010abb8f1"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def get_spectrum(state):\n",
    "    params = state.params['encoder']['stages_0']\n",
    "    lambda_bar = []\n",
    "    time_scales = []\n",
    "    for name, sequence_layer in params.items():\n",
    "        # read lambda parameters\n",
    "        Lambda_im = sequence_layer['S5SSM_0']['Lambda_im']\n",
    "        Lambda_re = sequence_layer['S5SSM_0']['Lambda_re']\n",
    "        \n",
    "        # read and compute delta and Lambda\n",
    "        delta = np.exp(sequence_layer['S5SSM_0']['log_step'][:, 0])\n",
    "        Lambda = Lambda_re + 1j * Lambda_im\n",
    "        \n",
    "        # compute lambda_bar and time scales\n",
    "        lambda_bar.append(np.exp(Lambda * delta))\n",
    "        time_scales.append(1 / np.abs(Lambda) / delta)\n",
    "    return lambda_bar, time_scales\n",
    "spectrum, time_scales = get_spectrum(state)"
   ],
   "id": "5a1602ed4a962265",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "Plot the spectrum of the recurrent operator and the corresponding time scales upon initialization.",
   "id": "2d1377b48ebf3728"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_spectrum(spectrum):\n",
    "    fig, axes = plt.subplots(1, 6, figsize=(len(spectrum) * 4, 4))\n",
    "    # draw the unit circle\n",
    "    theta = np.linspace(0, 2 * np.pi, 100)  # 100 points from 0 to 2*pi\n",
    "    x = np.cos(theta)\n",
    "    y = np.sin(theta)\n",
    "    \n",
    "    # plot the spectrum\n",
    "    for i, (ax, layer) in enumerate(zip(axes, spectrum)):\n",
    "        ax.plot(x, y, 'r', linewidth=1)\n",
    "        ax.scatter(np.real(layer), np.imag(layer), marker='o', alpha=0.8)\n",
    "    \n",
    "        # format axis\n",
    "        ax.set_title(f'Layer {i}')\n",
    "        ax.set_aspect('equal', adjustable='box')\n",
    "        ax.set_xlim(-1.1, 1.1)\n",
    "        ax.set_ylim(-1.1, 1.1)\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "plot_spectrum(spectrum)"
   ],
   "id": "9ff987be826d7314",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def plot_time_scales(time_scales):\n",
    "    log_scales = np.log2(np.stack(time_scales).flatten())\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n",
    "    ax.hist(log_scales)\n",
    "    \n",
    "    # format axis\n",
    "    max_scale = np.max(np.ceil(log_scales))\n",
    "    min_scale = np.min(np.floor(log_scales))\n",
    "    ax.set_xlim((min_scale, max_scale))\n",
    "    xticks = np.arange(1 + max_scale - min_scale) + min_scale\n",
    "    ax.set_xticks(xticks, (2 ** xticks).astype(np.int32))\n",
    "    ax.set_title('Distribution of time scales')\n",
    "    ax.set_xlabel('Time scale')\n",
    "    ax.set_ylabel('Count')\n",
    "    plt.show()\n",
    "    \n",
    "plot_time_scales(time_scales)"
   ],
   "id": "c1779ae6f5f72b44",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Train the model\n",
    "For training, we implemented a trainer module that makes training as easy as possible. The trainer module hides some boilerplate code for training from the user and provides a simple interface to train the model. It loops through the data loader, computes the loss, and updates the model parameters. Therefore, we need to define training_step and validation_step functions that the loop calls upon the model. These are implemented already, and can be used here."
   ],
   "id": "b4970c69f459df1d"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "from event_ssm.train_utils import training_step, evaluation_step\n",
    "from event_ssm.trainer import TrainerModule\n",
    "\n",
    "# just-in-time compile the training and evaluation functions\n",
    "train_step = jax.jit(training_step)\n",
    "eval_step = jax.jit(evaluation_step)\n",
    "\n",
    "# initialize the trainer module\n",
    "num_devices = 1\n",
    "trainer = TrainerModule(\n",
    "    train_state=state,\n",
    "    training_step_fn=train_step,\n",
    "    evaluation_step_fn=eval_step,\n",
    "    world_size=num_devices,\n",
    "    config=cfg,\n",
    ")"
   ],
   "id": "61ab72052c47f47c",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "We are now ready to start the training loop. \n",
    "\n",
    "**Note:** JAX compiles your program just-in-time (JIT) to optimize performance. This means that the first iteration of the training loop will be slower than the following ones.  "
   ],
   "id": "d66a413bc0ac7d2b"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# generate random key for dropout\n",
    "key, dropout_key = jax.random.split(key)\n",
    "\n",
    "# train the model\n",
    "trainer.train_model(\n",
    "    train_loader=train_loader,\n",
    "    val_loader=val_loader,\n",
    "    test_loader=test_loader,\n",
    "    dropout_key=dropout_key\n",
    ")"
   ],
   "id": "9d5ab8aa623db697",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Inspect the trained model\n",
    "We now have a trained toy model on the SHD dataset.\n",
    "Let's look into the spectrum of the recurrent operator after training."
   ],
   "id": "a929b74f8ce235e5"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "spectrum, time_scales = get_spectrum(trainer.train_state)\n",
    "plot_spectrum(spectrum)\n",
    "plot_time_scales(time_scales)"
   ],
   "id": "3281a08743303429",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Assignment\n",
    "The function `apply_ssm` in `event_ssm/ssm.py` implements the recurrent operator with an associative scan. On highly parallel GPUs, this can speed up training on very long sequences. \n",
    "On CPUs however, the overhead of the scan operation can slow down training. \n",
    "Your task is to implement a CPU-friendly version of the recurrent operator in `event_ssm/ssm.py` and compare the training time with the original implementation.\n",
    "We suggest to implement a step-by-step recurrence with [`jax.lax.scan`](https://jax.readthedocs.io/en/latest/_autosummary/jax.lax.scan.html) instead of the currenlty used [`jax.lax.associative_scan`](https://jax.readthedocs.io/en/latest/_autosummary/jax.lax.associative_scan.html) for this purpose."
   ],
   "id": "28ad17b8c7b61230"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}