how-to-train-your-perceptron-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "052ac8d4-9be9-4046-9a7f-7dc86881ff48",
   "metadata": {},
   "source": [
    "# How to Train your Perceptron-1\n",
    "**An opinionated recipe to get the most performance on your first few tries**\n",
    "\n",
    "For a more traditional recipe check out [A Recipe for Training Neural Networks](http://karpathy.github.io/2019/04/25/recipe/) by [Andrej Karpathy](https://karpathy.ai/)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e60995c-36c8-43e0-87bf-8cc3ee1785aa",
   "metadata": {},
   "source": [
    "### A perceptron (a.k.a. artificial neuron) is the building block of a traditional neural network (a.k.a Multilayer Perceptron; MLP)\n",
    "\n",
    "\n",
    "<img src=\"data/img/neural_net.png\">\n",
    "\n",
    "([image source](https://pubmed.ncbi.nlm.nih.gov/28087243/))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "559e00f9-b535-4061-aedd-5166dd2ddf47",
   "metadata": {},
   "source": [
    "In astronomy, we encounter mainly these three data types:\n",
    "- **Catalogs (i.e. tabulated numerical data)** &rarr; This Notebook\n",
    "- Images (i.e. data with spatial dependence/ evolution)\n",
    "- Time Series (i.e. data with temporal dependence/evolution)\n",
    "\n",
    "and we mainly encounter two types of tasks:\n",
    "- Regression (i.e. predict a continuous value)\n",
    "- Classification (i.e. predict a discrete value)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "680e59f1-7096-4e62-8cd8-8cb3bebf86ad",
   "metadata": {},
   "source": [
    "**Basic Training algorithm:**\n",
    "- Define a model (a sequence of matrix multiplication)\n",
    "- Start with an initial guess for all the parameters\n",
    "- Find the derivative of loss function with respect to parameters\n",
    "- modify the parameters by some step to reduce the derivate value.\n",
    "- take multiple steps until derivative is zero (ideally) or loss minima is reached.\n",
    "\n",
    "<img src=\"data/img/gradient_descent_animation.gif\">\n",
    "\n",
    "([image source](https://towardsdatascience.com/gradient-descent-animation-1-simple-linear-regression-e49315b24672#:~:text=The%20Gradient%20Descent%20method%20is,and%20costs%20during%20gradient%20descent.))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "673beee1-7e85-4a76-981e-39ce3d75a482",
   "metadata": {},
   "source": [
    "### Your first neural Network:\n",
    "Which Python framework to use?\n",
    "\n",
    "There are two main options, [Tensorflow](https://www.tensorflow.org/) and [Pytorch](https://pytorch.org/). Both give you the exact same functionalities.\n",
    "\n",
    "- In a hurry? &rarr; Use Tensorflow\n",
    "- Long research project? &rarr; Use Pytorch\n",
    "\n",
    "Today we are in a hurry!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3638af0d-67e4-4e47-96c9-c47bd9ef9aa9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import tensorflow as tf\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from tensorflow.keras import layers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb7445a1-e9d7-45d6-9408-35a38fc92962",
   "metadata": {},
   "source": [
    "We will take the case of photometric redshift estimation using magnitudes and colors as our example. This is a regression task. We will deal with classification exclusively in a future tutorial but will spread some tidbits for classification along the way. Let's read in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6d9f621-87c8-4a39-b0c9-14aa2550c443",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_path = Path(\"./data/photoz_catalogues/Teddy\")\n",
    "teddy_A = pd.read_csv(\n",
    "    base_path / \"teddy_A\",\n",
    "    delim_whitespace=True,\n",
    "    comment=\"#\",\n",
    "    names=[\n",
    "        \"id\",\n",
    "        \"mag_r\",\n",
    "        \"u-g\",\n",
    "        \"g-r\",\n",
    "        \"r-i\",\n",
    "        \"i-z\",\n",
    "        \"z_spec\",\n",
    "        \"feat1\",\n",
    "        \"feat2\",\n",
    "        \"feat3\",\n",
    "        \"feat4\",\n",
    "        \"feat5\",\n",
    "    ],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbc038da-4472-45c4-8571-4f2b2c3bc145",
   "metadata": {},
   "source": [
    "### Select features carefully and ALWAYS scale the data\n",
    "\n",
    "**NOTE:**\n",
    "Feature = predictors = independent variable = input | Target = output = dependent variable\n",
    "\n",
    "Neural Networks are not the best at automatic feature selection, so better to start with a small number of physically motivated features which are (hopefully) independent. \n",
    "\n",
    "Feature scaling takes away the effect of Units and does not let a feature dominate. There are many ways to scale, subtracting each feature's mean and dividing by their standard deviation aka Standard Scaling is a good start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c00c6552-6372-446d-957e-77f889eb5cd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = teddy_A[[\"mag_r\", \"u-g\", \"g-r\", \"r-i\", \"i-z\",]]\n",
    "Y = teddy_A[\"z_spec\"].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75ecf0a1-b49b-48b2-9c49-c9f522539670",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "X = scaler.fit_transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f7b11f0-a9ee-4740-a963-9fe2bbf54e32",
   "metadata": {},
   "source": [
    "## For regression: If the target variable is bounded scale it to an unbounded space using a sigmoidal function\n",
    "\n",
    "Bounded= values are either floored or capped or both\n",
    "\n",
    "\n",
    "One such function that works well is:\n",
    "$$h(x) = \\log \\left( \\dfrac{x-x_{\\mathrm{min}}}{x_{\\mathrm{max}}-x}\\right)$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63378f6b-26ee-4f34-a8ed-c3d24f3e4a9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistic_transform(x, xmin, xmax):\n",
    "\n",
    "    return np.log((x - xmin) / (xmax - x))\n",
    "\n",
    "\n",
    "def inverse_logistic_transform(x, xmin, xmax):\n",
    "\n",
    "    return (xmax * np.exp(x) + xmin) / (1 + np.exp(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53269757-793f-4bc6-bbfd-9c8b48eec7f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = logistic_transform(Y, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e055f16b-e738-40a1-b734-46e3073b24c4",
   "metadata": {},
   "source": [
    "### Set a random seed\n",
    "Much of the following will be probabilistic in nature. For the sake of reproducibility during development se a random seed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e975fe4a-810c-4ac6-b89b-2d3f6fab7a79",
   "metadata": {},
   "outputs": [],
   "source": [
    "SEED = 42\n",
    "tf.random.set_seed(SEED)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a540691b-df80-4491-832e-27ae7c6ac1ae",
   "metadata": {},
   "source": [
    "### Randomly divide available data into 3 sets: Training, validation, test\n",
    "\n",
    "There is a fourth data set which is not available to us at this moment. This will be the one used to evaluate the method when everything has been finalised. If the data set is large, 80:10:10 is a good split."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b303d92-0671-4973-8647-4426c98398d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_gal = len(X)\n",
    "frac_train = 0.8\n",
    "frac_val = 0.1\n",
    "\n",
    "rng = np.random.default_rng(SEED)\n",
    "indices = rng.permutation(n_gal)\n",
    "\n",
    "ind_split_train = int(np.ceil(frac_train * n_gal))\n",
    "ind_split_val = ind_split_train + int(np.ceil(frac_val * n_gal))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5637cd4e-f5ef-42ed-8652-5fca8123d99c",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = X[indices[:ind_split_train]]\n",
    "x_val = X[indices[ind_split_train:ind_split_val]]\n",
    "x_test = X[indices[ind_split_val:]]\n",
    "\n",
    "y_train = Y[indices[:ind_split_train]]\n",
    "y_val = Y[indices[ind_split_train:ind_split_val]]\n",
    "y_test = Y[indices[ind_split_val:]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce0e79a4-f61f-4638-bdb4-189231338b50",
   "metadata": {},
   "source": [
    "Choosing MLP architecture:\n",
    "- Start small and increase model size untill performance saturates for a fixed amount of data\n",
    "- Wider networks are good at memorization whereas deep networks are good at generalization\n",
    "- Generally a wide first- narrow later structure works well regression\n",
    "\n",
    "\n",
    "Choosing activation function: (See [this](https://github.com/christianversloot/machine-learning-articles/blob/main/how-to-use-prelu-with-keras.md) for reasoning)\n",
    "- Always use the \"leaky reLu\" family for hidden layers\n",
    "- Use linear activation for the output layer for regression and sigmoid/softmax for classification\n",
    "\n",
    "Choosing an initializer for weights:\n",
    "- For reLu family of activations use He-normal\n",
    "- For Sigmoid/Tanh family of activations use Glorot-Uniform\n",
    "\n",
    "See [this](https://www.deeplearning.ai/ai-notes/initialization/) to see why initialization is important. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2848621-66d0-42bb-87f1-8bb7fd398636",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_input = layers.Input(shape=(5,))\n",
    "h1 = layers.Dense(100, kernel_initializer=\"he_normal\")(model_input)\n",
    "o1 = layers.PReLU()(h1)\n",
    "h2 = layers.Dense(50, kernel_initializer=\"he_normal\")(o1)\n",
    "o2 = layers.PReLU()(h2)\n",
    "output = layers.Dense(1, activation=\"linear\")(o2)\n",
    "model = tf.keras.models.Model(inputs=model_input, outputs=output)\n",
    "# summarize layers\n",
    "print(model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "529b964c-cebc-4491-ac1c-f33e81e56ed1",
   "metadata": {},
   "source": [
    "**Choosing an optimizer:**\n",
    "- Just use Adam with default learning rate\n",
    "\n",
    "**Choosing a loss**\n",
    "- Mean squared error (MSE) for regression and cross entropy for classification. They maximize $p(\\mathrm{data}|\\mathrm{parameters})$ under the assumption of a gaussian data generating process and are good starting points.\n",
    "- Classification accuracy is a good metric for classification. MSE is a good metric for regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25d09d44-1d2b-4279-b674-2c26f4726436",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.compile(\n",
    "    optimizer=tf.keras.optimizers.Adam(learning_rate=1e-3),\n",
    "    loss=tf.keras.losses.MeanSquaredError(),\n",
    "    metrics=tf.keras.losses.MeanSquaredError(),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdf841f3-8d70-4718-a398-a27389dfcf8d",
   "metadata": {},
   "source": [
    "**Fitting the model**\n",
    "- Batch size: The general norm is to maximize it until you run out of memory. This is a good recipe for image data but not for catalogs, where often the whole data set can fit in memory. In that case the network will see only one batch and overfit to the training data.\n",
    "- Shuffle the data during each epoch. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31d2c653-fa06-4b73-8f92-7d183c699952",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_epochs = 100\n",
    "batch_size = 500\n",
    "\n",
    "history = model.fit(\n",
    "    x=x_train,\n",
    "    y=y_train,\n",
    "    batch_size=batch_size,\n",
    "    epochs=n_epochs,\n",
    "    verbose=1,\n",
    "    callbacks=None,\n",
    "    validation_data=(x_val, y_val),\n",
    "    shuffle=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc5dbcbe-5aed-4b48-b7ad-4106dff7c6d6",
   "metadata": {},
   "source": [
    "**How many epochs do we let it train for?**\n",
    "\n",
    "\n",
    "Plot the train and validation loss as a function of epoch\n",
    "- If everything is right, training loss should decrease with increasing epoch and validation loss decreases reaches a minima and then increases.\n",
    "- We select the model with the minimum validation loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4e5a96e-3d94-4212-ae89-84419b9786b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "epoch = np.arange(1, n_epochs+1)\n",
    "\n",
    "plt.plot(epoch, history.history[\"loss\"], label=\"Training loss\")\n",
    "plt.plot(epoch, history.history[\"val_loss\"], label=\"Validation loss\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75d5a904-3c11-4f6b-b6ec-b25f73b234f4",
   "metadata": {},
   "source": [
    "**If performance above looks good, we are done! Else:**\n",
    "- Learning rate decay\n",
    "- Change the model architecture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b584c4b6-bad5-41bc-9617-741d9b7b2762",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tf2",
   "language": "python",
   "name": "tf2"
  },
  "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.9.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}