08_08_NumPy_CrossEntropy_Softmax_KL_Mean_Std.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5859fd3b",
   "metadata": {},
   "source": [
    "# Learning Objectives\n",
    "\n",
    "\n",
    "Replacing inefficient Pyth on loops, large trip count loops that perform low level operations, will negativley impact numeric computation performance.\n",
    "\n",
    "Replacing such loops has the following benefits, which are related to this modules learning objectives.\n",
    "\n",
    "| PyTorch UFUNCs | NumPy UFUNCs | Description |\n",
    "| --- | --- | --- |\n",
    "| torch.sigmoid(x)|1 / (1 + np.exp(-x))| **Sigmoid Function** |\n",
    "| torch.nn.Softmax() | np.exp(x)/np.exp(x).sum()| **Softmax function** |\n",
    "| torch.nn.CrossEntropyLoss() | 1.0 / (1.0 + np.exp(-x))| **Cross Entropy Loss** |\n",
    "| torch.nn.KLDivLoss() | np.sum(p * np.log(p/q) | **Kullback-Leibler Divergence (KL)** |\n",
    "\n",
    "Code will be:\n",
    "    \n",
    "- More **readable**\n",
    "- More **maintainable** for the developer\n",
    "- Maintained by 3rd parties who invest THEIR time perfecting the algorithm\n",
    "- **Faster** on existing hardware\n",
    "- Faster on future enhancements to hardware platforms and 3rd party library developments\n",
    "\n",
    "### Learning Objectives\n",
    "At the end of this module you will be able to:\n",
    "- Apply NumPy vectorized libraies for inefficient loopy code\n",
    "- Describe the benefits of using NumPy as an alternative to your \"roll your code\"\n",
    "\n",
    "**To run the lab**: These step could be run on a laptop - NOT REQUIRED for DevCloud!\n",
    "**Laptop Requirements:**\n",
    "\n",
    "```bash\n",
    "conda config --add channels intel\n",
    "conda install numpy\n",
    "conda install scipy\n",
    "conda install update pandas\n",
    "```\n",
    "\n",
    "## Python loops are bad for performance\n",
    "\n",
    "**Python is great!** Its a great language for AI. There are many, many advantages in using Python especially for data science.\n",
    "\n",
    "- Easy to program (don’t worry about data types and fussy syntax at least relative to C/C++ and other languages\n",
    "- FAST for developing code!\n",
    "- Leverages huge array of libraries to conquer any domain\n",
    "- Lots of quick answers to common issues in Stack Exchange\n",
    "\n",
    "**Python, however, is slow for Massively repeating small tasks** - such as found in loops! \n",
    "-Python loops are SLOW\n",
    "-Compared to C, C++, Fortran and other typed languages\n",
    "-Python is forced to look up every occurrence and type of variable in a loop to determine what operations it can perform on that data type\n",
    "-It cannot usually take advantage of advances in hardware in terms of vector width increases, multiple cores, new instructions from a new HW instruction set, new AI accelerators, effective cache memory layout, and more\n",
    "\n",
    "**BUT: Python has library remedies to these ills!**\n",
    "\n",
    "Importing key libraries shift the burden of computation to highly efficient code.\n",
    "\n",
    "**NumPy**, for example, through its focus on elementwise efficient operations, gives indirect access to the efficiencies afforded in \"C\"\n",
    "\n",
    "libraries included in oneAPI and NumPy, SciPy, Scikit-learn all powered by Intel(r) oneAPI give access to modern advancements in hardware level: access to better cache and memory usage, access to low level vector instructions, and more.\n",
    "\n",
    "By leveraging packages such as these powered by oneAPI AND keeping libraries up to date, more capability is added to your underlying frameworks so that moving code, especially in a cloud world, can give you ready access to hardware acceleration, in many cases, without having to modify code this vectorized code\n",
    "Routines are written in C (based on Cython framework)\n",
    "\n",
    "NumPy arrays are densely packed arrays of homogeneous type. \n",
    "\n",
    "Python lists, by contrast, are arrays of pointers to objects, even when all of them are of the same type. So, you get the benefits of not having to check data types, and you also get locality of reference. Also, many NumPy operations are implemented in C, avoiding the general cost of loops in Python, pointer indirection and per-element dynamic type checking. The speed boost depends on which operations you’re performing.\n",
    "\n",
    "Goal of this module: **Search and destroy (replace) loops**\n",
    "\n",
    "- Avoid loops if you can - find an alternative if possible. \n",
    "- Sometimes it cannot be done - true data dependencies may limit our options. But many, many time there are alternatives.\n",
    "\n",
    "**The problem**\n",
    "\n",
    "- Loops isolate your code from hardware and software advances that update frequently.\n",
    "- They prevent you from effectively using key underlying resources - it is a waste.\n",
    "- They consume your time!\n",
    "- They can waster energy!\n",
    "\n",
    "\n",
    "## Reference:\n",
    "\n",
    "Video: Losing your Loops Fast Numerical Computing with NumPy by Jake VanderPlas .\n",
    "\n",
    "Book: Python Data Science Handbook by Jake VanderPlas.\n",
    "\n",
    "Book: Elegant SciPy: The Art of Scientific Python by by Juan Nunez-Iglesias, Stéfan van der Walt, Harriet Dashnow\n",
    "\n",
    "Article: The Ultimate NumPy Tutorial for Data Science Beginners : by Aniruddha April 28, 2020 at www.analyticsvidhya.com\n",
    "\n",
    "Academic Lecture pdf: Vectorization by Aaron Birkland Cornell CAC\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de622d6f",
   "metadata": {},
   "source": [
    "# Exercise: Compute Mean & Std of array\n",
    "\n",
    "Below is an example of a loop based way to compute the mean and standard deviation for a list or vector of values\n",
    "\n",
    "```python\n",
    "for i in range (len(a)):\n",
    "    S += a[i]\n",
    "mean = S/len(a)\n",
    "std = 0\n",
    "for i in range (len(a)):\n",
    "    d = a[i] - mean\n",
    "    std += d*d\n",
    "std = np.sqrt(std/len(a))\n",
    "print(\"mean\", mean)\n",
    "print(\"std\", std)\n",
    "```\n",
    "\n",
    "In the following exercise, replace this code with a more readbale and maintainable NumPy vectorized variant as follows:\n",
    "\n",
    "```python\n",
    "print(a.mean())\n",
    "print(a.std())\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7b83ebd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import time\n",
    "import math\n",
    "\n",
    "rng = np.random.default_rng(2021)\n",
    "# random.default_range is the recommended method for generated random's\n",
    "# see blog \"Stop using numpy.random.seed()\" for reasoning\n",
    "# https://towardsdatascience.com/stop-using-numpy-random-seed-581a9972805f\n",
    "\n",
    "a = rng.random((10_000_000,))\n",
    "a_torch = torch.from_numpy(a)\n",
    "t1 = time.time()\n",
    "timing = {}\n",
    "S = 0\n",
    "\n",
    "################################ code to replace in next cell ##############\n",
    "for i in range (len(a)):\n",
    "    S += a[i]\n",
    "mean = S/len(a)\n",
    "std = 0\n",
    "for i in range (len(a)):\n",
    "    d = a[i] - mean\n",
    "    std += d*d\n",
    "std = np.sqrt(std/len(a))\n",
    "print(\"mean\", mean)\n",
    "print(\"std\", std)\n",
    "############################################################################\n",
    "timing['loop'] = time.time() - t1\n",
    "\n",
    "\n",
    "print(timing)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c92aa67c-6e5b-4156-9dac-b4d6b945f1c1",
   "metadata": {},
   "source": [
    "# Show Intel oneMKL-DNN under the hood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8168c553-86f2-4341-864b-bc146b29f859",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(torch.__config__.parallel_info())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b12d8d2e",
   "metadata": {},
   "source": [
    "## Excercise: use NumPy mean and std\n",
    "This cell will error - fix the error\n",
    "\n",
    "Hint:\n",
    "\n",
    "```python\n",
    "print(a.mean())\n",
    "print(a.std())\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d329b1f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = time.time()\n",
    "##### insert NumPy code here ###############\n",
    "#print(np.xxx())\n",
    "print(a.mean())\n",
    "print(a.std())  \n",
    "############################################\n",
    "\n",
    "timing['numpy'] = time.time() - t1\n",
    "print(timing)\n",
    "print(f\"NumPy Acceleration {timing['loop']/timing['numpy']:4.1f} X\")\n",
    "\n",
    "t1 = time.time()\n",
    "##### insert NumPy code here ###############\n",
    "#print(np.xxx())\n",
    "print(a_torch.mean())\n",
    "print(a_torch.std())  \n",
    "############################################\n",
    "\n",
    "timing['numpy'] = time.time() - t1\n",
    "print(timing)\n",
    "print(f\"PyTorch Acceleration {timing['loop']/timing['numpy']:4.1f} X\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2294b5a-db88-4c20-b5ee-32edc80bab04",
   "metadata": {},
   "source": [
    "```python\n",
    "rng = np.random.default_rng(2021)\n",
    "a = rng.random((10_000_000,))\n",
    "a = rng.random((10_000_000,))\n",
    "a_torch = torch.from_numpy(a)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c0f1a4c-31cf-49ad-a43d-be99d2557178",
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = time.time()\n",
    "##### insert NumPy code here ###############\n",
    "#print(np.xxx())\n",
    "a.mean()\n",
    "a.std()\n",
    "############################################\n",
    "\n",
    "timing['numpy'] = time.time() - t1\n",
    "print(f\"NumPy Acceleration   {timing['loop']/timing['numpy']:4.1f} X\")\n",
    "\n",
    "t1 = time.time()\n",
    "##### insert NumPy code here ###############\n",
    "#print(np.xxx())\n",
    "a_torch.mean()\n",
    "a_torch.std() \n",
    "############################################\n",
    "\n",
    "timing['numpy'] = time.time() - t1\n",
    "print(f\"PyTorch Acceleration {timing['loop']/timing['numpy']:4.1f} X\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db1f9617",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.title(\"Measure acceleration of looping versus Numpy Mean & STD [Lower is better]\",fontsize=12)\n",
    "plt.ylabel(\"Time in seconds\",fontsize=12)\n",
    "plt.yscale('log')\n",
    "plt.xlabel(\"Various types of operations\",fontsize=14)\n",
    "plt.grid(True)\n",
    "plt.bar(x = list(timing.keys()), height= list(timing.values()), align='center',tick_label=list(timing.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e25d7aa",
   "metadata": {},
   "source": [
    "# Cross Entropy\n",
    "\n",
    "Cross entropy calculates are done all the time in machine learning. consider  the loopy stream of consciousness variant of computing cross entropy below, as seen in the next cell. Try your hand at removing the loop and using the NumPy Ufunc, aggregation or other NumPy construct to make this code more readable and faster.\n",
    "\n",
    "###  Logistic Regression\n",
    "\n",
    "probability of y=1\n",
    "\n",
    "# $ \\hat{y} = q_{(y=1)} = g(w,x) = \\frac{1}{1 + e^{-wx}} $ \n",
    "\n",
    "probability of y = 0\n",
    "\n",
    "\n",
    "# $ q_{(y=1)} = 1-\\hat{y}$\n",
    "\n",
    "w = weights\n",
    "\n",
    "x = input vector\n",
    "\n",
    "\n",
    "## cross-entropy to get a measure of dissimilarity between p and q\n",
    "\n",
    "# $ H(p,q) =  - \\Sigma p_i log(q_i) = -y log \\hat{y} - (1-y) log(1-\\hat{y}) $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10f483f9",
   "metadata": {},
   "source": [
    "# Slow python method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed0f384d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "from math import log2\n",
    "import numpy as np\n",
    "import time\n",
    "# Sigmoid function\n",
    "\n",
    "timing = {}\n",
    "\n",
    "######################### This is the targted function for this exercise ######\n",
    "def sigmoid_slow(z):\n",
    "    sig = [1.0/(1.0+math.exp(-1*xi)) for xi in z]\n",
    "    return sig\n",
    "###############################################################################\n",
    "\n",
    "\n",
    "# yHat represents the predicted value / probability value calculated as output of hypothesis / sigmoid function\n",
    " \n",
    "# y represents label\n",
    "\n",
    "# cross_entropy list comprehension method\n",
    "\n",
    "def cross_entropy_loss_slow(yHat, y):\n",
    "    if y == 1:\n",
    "        return [-1*math.log(yi) for yi in yHat]\n",
    "    else:\n",
    "        return [-1*math.log(1.0 - yi) for yi in yHat]\n",
    "   \n",
    "x = [xi/100_000. for xi in range(-1_000_000, 1_000_000)]  # num between -10 and 10 step .00001\n",
    "p = x\n",
    "q = []\n",
    "start = time.time()\n",
    "sig_x = sigmoid_slow(x)\n",
    "cost_1 = cross_entropy_loss_slow(sig_x, 1)\n",
    "cost_0 = cross_entropy_loss_slow(sig_x, 0)\n",
    "timing['list comprehension'] = time.time() - start\n",
    "print(f\"time timing: {timing['list comprehension']:5.3f} sec\")\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "plt.plot(sig_x, cost_1, label='$ \\hat{y} $ if y=1')\n",
    "plt.plot(sig_x, cost_0, label='$ (1 - \\hat{y} )$ if y=0')\n",
    "plt.xlabel('sigmoid(x)')\n",
    "plt.ylabel('$ \\hat{y}$')\n",
    "plt.legend(loc='best')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f79978c9",
   "metadata": {},
   "source": [
    "## Hint:\n",
    "\n",
    "```python\n",
    "def sigmoid(x):\n",
    "    return 1.0 / (1.0 + np.exp(-x))\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f934884",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import torch\n",
    "# Sigmoid function NumPy\n",
    "def sigmoid(z):\n",
    "##### insert improved NumPy code here\n",
    "    return 1.0 / (1.0 + np.exp(-z))\n",
    "#####################################\n",
    "\n",
    "\n",
    "# yHat represents the predicted value / probability value calculated as output of hypothesis / sigmoid function\n",
    " \n",
    "# y represents label\n",
    "\n",
    "# cross_entropy NumPy method\n",
    "\n",
    "def cross_entropy_loss(yHat, y):\n",
    "    if y == 1:\n",
    "        return -np.log(yHat)\n",
    "    else:\n",
    "        return -np.log(1 - yHat)\n",
    "    \n",
    "x = np.arange(-10, 10, 0.00001)\n",
    "start = time.time()\n",
    "sig_x = sigmoid(x)\n",
    "cost_1 = cross_entropy_loss(sig_x, 1)\n",
    "cost_0 = cross_entropy_loss(sig_x, 0)\n",
    "\n",
    "timing['numpy'] = time.time() - start\n",
    "print(f\"time elapsed: {timing['numpy']:5.3f} sec\")\n",
    "\n",
    "ratio = timing['list comprehension']/timing['numpy']\n",
    "print(f'Acceleration: {ratio:5.1f}X')\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "plt.plot(sig_x, cost_1, label='$ \\hat{y} $ if y=1')\n",
    "plt.plot(sig_x, cost_0, label='$ (1 - \\hat{y} )$ if y=0')\n",
    "plt.xlabel('sigmoid(x)')\n",
    "plt.ylabel('$ \\hat{y}$')\n",
    "plt.legend(loc='best')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c56cea94-5fae-4f09-9865-3dcfbffe33f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.title(\"Measure acceleration of looping versus Numpy Cross Entropy Loss [Lower is better]\",fontsize=12)\n",
    "plt.ylabel(\"Time in seconds\",fontsize=12)\n",
    "plt.yscale('log')\n",
    "plt.xlabel(\"Various types of operations\",fontsize=14)\n",
    "plt.grid(True)\n",
    "plt.bar(x = list(timing.keys()), height= list(timing.values()), align='center',tick_label=list(timing.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "092fee78-fa99-42f9-a7eb-8e92ba264c4f",
   "metadata": {},
   "source": [
    "# Implement a PyTorch CrossEntropyLoss\n",
    "\n",
    "The cell below will error - replace the xxx with valid code!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9db87ab7-4b47-4288-a9e7-01ffaeed0bdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import torch\n",
    "\n",
    "timing = {}\n",
    "\n",
    "start = time.time()\n",
    "x = np.arange(-10, 10, 0.00001)\n",
    "numpy_sig_x = sigmoid(x)\n",
    "timing['sigmoidNumPy'] = time.time() - start\n",
    "print(numpy_sig_x)\n",
    "del x\n",
    "\n",
    "x = torch.arange(-10, 10, 0.00001)\n",
    "start = time.time()\n",
    "torch_sig_x = torch.sigmoid(x)\n",
    "timing['sigmoidPyTorch'] = time.time() - start\n",
    "print(torch_sig_x.cpu().detach().numpy())\n",
    "\n",
    "ratio = timing['sigmoidNumPy']/timing['sigmoidPyTorch']\n",
    "print(f\"\\nAcceleration PyTorch to NumPy: {ratio:4.1f} X\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8065bc7-fabd-492f-becc-923d09ff0629",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import torch\n",
    "\n",
    "# yHat represents the predicted value / probability value calculated as output of hypothesis / sigmoid function\n",
    " \n",
    "# y represents label\n",
    "\n",
    "# cross_entropy NumPy method\n",
    "\n",
    "\n",
    "x = torch.arange(-10, 10, 0.00001,  requires_grad=True)\n",
    "\n",
    "##############  Use PyTorch to compute CrossEntropyLoss  ###############\n",
    "\n",
    "xxx \n",
    "\n",
    "########################################################################\n",
    "\n",
    "timing['torch'] = time.time() - start\n",
    "print(f\"time elapsed: {timing['torch']:5.3f} sec\")\n",
    "\n",
    "ratio = timing['list comprehension']/timing['torch']\n",
    "print(f'Acceleration: {ratio:5.1f}X')\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "plt.plot(sig_x, cost_1, label='$ \\hat{y} $ if y=1')\n",
    "plt.plot(sig_x, cost_0, label='$ (1 - \\hat{y} )$ if y=0')\n",
    "plt.xlabel('sigmoid(x)')\n",
    "plt.ylabel('$ \\hat{y}$')\n",
    "plt.legend(loc='best')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0a8db45",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.title(\"Measure acceleration of looping versus Numpy Cross Entropy Loss [Lower is better]\",fontsize=12)\n",
    "plt.ylabel(\"Time in seconds\",fontsize=12)\n",
    "plt.yscale('log')\n",
    "plt.xlabel(\"Various types of operations\",fontsize=14)\n",
    "plt.grid(True)\n",
    "plt.bar(x = list(timing.keys()), height= list(timing.values()), align='center',tick_label=list(timing.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e08e9e72",
   "metadata": {},
   "source": [
    "# Softmax Loop\n",
    "\n",
    "Another Algorithm that's used all the time in machine learning is Softmax.\n",
    "\n",
    "The softmax function, or normalized exponential function, converts a vector of K real numbers into a probability distribution of K possible outcomes.\n",
    "\n",
    "### Below is slower python for loop method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a00b0c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import time\n",
    "import math\n",
    "from sys import getsizeof\n",
    "\n",
    "timing = {}\n",
    "np.random.seed(seed=42)\n",
    "BIG = 10_000_000\n",
    "b = list(np.random.rand(BIG))\n",
    "def softmax_slow(x):\n",
    "    denominator = 0.0\n",
    "    for xi in x:\n",
    "        denominator += math.exp(xi)\n",
    "    return [math.exp(xi)/denominator for i, xi in enumerate(x)]\n",
    "start = time.time()\n",
    "softmax_slow(b)\n",
    "timing['softmax_loop'] = time.time() - start\n",
    "print(f\"time elapsed: {timing['softmax_loop']:5.3f} sec\")\n",
    "#print(f'memory: {getsizeof(b):,}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bbe6945",
   "metadata": {},
   "source": [
    "# Softmax NumPy & PyTorch\n",
    "\n",
    "More Readable/Maintainable/Faster numpy method:\n",
    "\n",
    "```python\n",
    "def softmax(x):\n",
    "    return(np.exp(x)/np.exp(x).sum())  # one line of code, no loop indices \n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "894fd136",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=42)\n",
    "b = np.random.rand(BIG)\n",
    "b_torch = torch.from_numpy(b)\n",
    "\n",
    "def softmax_numpy(x):\n",
    "    ########### insert solu\"tion here\n",
    "    return(np.exp(x)/np.exp(x).sum())  # one line of code, no loop indices \n",
    "    ################################\n",
    "    \n",
    "def softmax_torch(x):\n",
    "    #doIt = torch.nn.Softmax(dim=0)\n",
    "    #return doIt(x)\n",
    "    return(torch.exp(x)/torch.exp(x).sum())\n",
    "\n",
    "\n",
    "start = time.time()\n",
    "softmax_numpy(b)\n",
    "timing['softmax_numpy'] = time.time() - start\n",
    "print(f\"time elapsed: {timing['softmax_numpy']:5.3f} sec\")\n",
    "ratio = timing['softmax_loop'] / timing['softmax_numpy'] \n",
    "print(f'NumPy Acceleration: {ratio:5.4g}X')\n",
    "print(softmax_numpy(b))\n",
    "\n",
    "    \n",
    "start = time.time()\n",
    "softmax_torch(b_torch)\n",
    "timing['softmax_pytorch'] = time.time() - start\n",
    "print(f\"time elapsed: {timing['softmax_pytorch']:5.3f} sec\")\n",
    "ratio = timing['softmax_loop'] / timing['softmax_pytorch'] \n",
    "print(f'PyTorch Acceleration: {ratio:5.4g}X')\n",
    "print(softmax_torch(b_torch))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ea116cf-975f-43c8-81d0-f7a5055874ff",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c36f9746",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.title(\"Measure acceleration of looping versus Numpy Softmax [Lower is better]\",fontsize=12)\n",
    "plt.ylabel(\"Time in seconds\",fontsize=12)\n",
    "plt.yscale('log')\n",
    "plt.xlabel(\"Various types of operations\",fontsize=14)\n",
    "plt.grid(True)\n",
    "plt.bar(x = list(timing.keys()), height= list(timing.values()), align='center',tick_label=list(timing.keys()))"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAAyCAYAAAC+qXUzAAAREElEQVR4nO2deZxWVRnHvy8zgwwgMCDgQioqDm7gCrng4IILkLtkmJWm5kJhmkuWJYqSpqFmmlupqbmgolZgarhgLlhuIGqGS7nmvqYo0x+/e7r3Pe+56/u+4gzn+/nweZl7zz13O+c5z3mWc8Hj8Xg8Ho/H4/F0LPYDNspRvitwKtCjlhfRWMvKPJ5liPFAX+DymP3DgHNi9n0EXAX8ztr+NWBjx/YkPgFuBm4EdgY+iym3P/AQ8HiOuj0eTw4moE7WPaVcL+AJYCHQEGzrAuwKLAYOipRtARYAfQpcTwm4EDg8oUx3YBawXoH6PR5PCoOBZ4AVM5TtjgTAeY59DwMvoE4NcDowrYrragX+iQRPHKsC84FVqjiPx+OxaAQeBDbPWH4M0A5MtLY3AR8irQMkHF4DNqvi2kpIMGydUm5/4LYqzuPxeCzOCP5l5RQkGFawtk8Itn89+HsDZCtw2fxWQRrHDGB7ZNuYAezkKHsF8JMM13UvcGCGch6PJ4VW4A1g+RzHzAXmodG8C7AWcADwIuXThgnAk47jVwKuAUYiu8QHwNXAdkjjsG0cU4L9aWwYHD84roD3Sng82TgGuAh4L2P5ZmAE8CpwJ+rUC5G6Px7ZGAyDgDcddRwLTAZeQVOE7sDxSEjcirSMKO9SqZ24eAT4C3Jzfi3T3Xg8ngpWBj4GvpTjmB3QdOHbGcoejYRHEj8CniU0Vro4GPh7pquDUcASYIhrp9cYOjYjUWBLO8kNJiuP4B65lnWOQEbHf+U4pi34nZWh7NtAt5Qyo4A70LuOoxcSYFmYCzwNTEJaSRleMHRsjgL2Dv6/ELgupbwRIC1oFByE5pvLBft/EdTpCVkeOASYmvO4NtTxXs5Q9jWgt2P7UWhUvxDYhjBOoTua1uxrle+T8XygtnAtMkIeQbLA8XQwWtAo1o4a0KgCdSyHRo2XgbeocWhtJ2Af9Hw3zFh+CLAbikC8AxhNujawOop3aLC2vwacBByK3s0ewfapwC6Oeq5FdoOstKF728TekRQM4fni8xYyHhlN4PdIWOThY+BcYG1kQf9WLS+wEzAChTA/mrH8V1Cuw1SkrmcRDM8hjc/OkfgGerefohiH0chmcQ0Kg45SQgNDmq0iyoNB/TtkKVwinK+WrL/rTTfgAuQaclGi2DWNIbz5nsDPCl5fbyRlV3PsGwlcnPO6asUJ6AW3AzdVUU8vYDbp95C3jXwZODvntexI9piBppx152EusuDXm8OAk6s4fiTwPPnMAyXgKTK0mTOQ8akduVfuQhJoEfBnJLHqyS+pjBID+W0fJ5Sedwb//opcNKeR3DD3QC4eUOLL9JzXtSuav08J6roe+BOVo/Nk4Oc5664FJWAOoXCYVEVdq5I8nZiOjGXtaBpzJ2onz6I59VZW+RWAx0gfNV2cQrLN45bgOu4qUHcWGlA/OKtO9UdpRm08T5yEoQRchrwSeZmB3mPqgLY7etjbWyfeC/lwv1Pg5FnYFKk2cTQC7wCXWtvXRsLiBwnH7k44PxtA9hfdHalss6k0Dl2ErPjReWEDin9fI2P9tWQl4D/o3X2MounqxVeD89g2jVOB95FwMVyKVOIi9EBxAH0TylxJPvU5Dxuj+/xeneq32Qw4scBxuwC/KXjOU9E9piWEcS4yhLhGjWlIo6iHgeoPJM9vN0M3cIhj33tIe4hjTyTYAAaSTa3tBtyDAlFc6tm6wfXsb20/FvhVhvrrwc6EWsOT1M+QeAFqI3ZjMoY6E247BHXsalT9X6CpUhy/Au6uov4k9kL3M65O9bvYARieo3w34tO7szAJ3eOw6EaX8XEUGrk/cOx7DKnPeS48C31Q7PcfE8qYaYw9OqyD7AbzU86Rd+5/Dpob74s0Epvng1/bonsX0lCWhq1hFnBJ8P9WNDWrNSUUhTcPhdVGGRH8Lgh+JyDL/OIqzncXobbnIs7Nlvf5u9KRjVr/Us66quHPZDd0AvyX6jQaE5sR1fIqBMMKSHLEqWb9gl+Xz7UaNkcjy+sJZdqQmm7HlB+KGmiSn9nEqkM2f+36KFf+Tyif3oUJPbVH5YeDfc6Iss+BycjCDdJmdk0oW4QBwFA0vYrSE60+dBNwX7BtSzTdsumOFiN5ACX0jEODwm1UGp7vQ22yV8I1Rd9pCRnyrkCuvluQgIqyOxI4c9Dc/Gr0ru20Z3POtxPO3dEx9zYwutFWkU3KZpwV1mgKj6OR8in08HqhDvsl9JI+Qn7cQajDL0YJG0/gjswagmLI4zptQ3Bt9wPbok6+NrKDbIKCP16IObYI+wW/tj0jimnAD1nbP0HaxFBkjHMxBvhxxmu5Gjg/Y1mQprcnGnWagN8iQVerUc+0kffRc++FRtv9UbuJhgCvTajBRDkPGbofRR33PDQ4zEYGtKjB8T/ontZF799FVDs4FsURbI3eRT8knJqQPWJzZDxeCxnVT0JToDGovUYxGkNnFgxG6ysTvLZgGI3U5nsdFZSAsegh/Rt1nuuQD/ZCJBg2QurrPOBIZKE+CHXsmehFuATDANTQ4hiOXtInSHN4F42KxyGBErecVRxp8Rum09+TUGbn4PcBx74PSTaY3QbcnnINhiIRaQtR0s90NPW7hPB6q2Wb4LcP6nz/RnP8S9AgYCghzcl+rysie5BRl1dBlvG3gRtQp7VJe57mGTUh490JhAlGb6DAn6lIMGyFgsEWBfufRgPT21QKctNZ7CmTi5FIa/oi8RbpuROmP3aNbrQFQxvqDP91VLA1CqM1Kvs05Ea8A41qjUjl2xi9DIBfB393Jdna+hnJnXvb4PdH5Jt/GZYQjipZOlpf1LDeiNm/PBKMs6nUGIxLNc09V+8Q1LOQIB+DRsceuO1GeSihd/E3suX9m2cR5RU0kIA63rCgrg9xtxFTR3PMOaJCvhVFctoj/EcourAFCeTTkTD4B7AmEqSvOeo277CJ9ByEA0lIY15KPEa6YDD2n7LBMioY+qGX5Ho5jcg6/DSVQRhd0JzxfOT6cHWmtE7wJskdqQ29uCJCAXSNdkBOEo8jI2xP5PHoCfRHvnqQuroYtxuuhATLuwn1b0CY45DGXGSQKsJxyCA4nuqFAkizayX7EmRvEd+hQRpqF8I4hObg7+i1mtyOd2LqiLYtY6Oy7T79kKB/Dxmp/4Cmcr2RljM2pu4lwW8LyRotlK/f2JEwWtGS6MaoYDCjsh0s0h+5p1ZAGoKdA96MVrwdjNR6m3bSVf0n0dy0RKUQaUINKEuWWhJ5RuhLkQFrGyQMjkSjzkxCd9zuaP7roi9ajCOOF8nue38mYzmb3ki93xHZgmpBXBtxYVymA63tbWgKOhqFDy8g7HTTkTcoavDthgaeuOfZQPhuX0Ea3IjIfuNFmUUYWrwIR0ahAzOF6E2+zMqOhHEklA0cjWh6cDBq6CDpORq9jDWRseYJ9EBd6tYENH+9Fs1rXUkc9gi9GqG7D2R5XpNKwTApOG9P1MCmBP/KpFtG8hwzDy2IcRpqEOOQAJyLbCobokboYmUkhZO0mzepb5htV+Qd+AluG0heBiFhuGfw906ofdyYctx9yAgbZQgyhPZD9oZm1PlHIW3T9gKtgzqoS7gdgQLjBiBNdxoygs4GfoqWdj8daSEmGnQBarMbEQ5YbwBnEnpTDKa9J9k3OjpGMJRp+o1oRG5Ho+FMJDleD/7dgCy6LpsDqCNfjjSFg5F9Ygbx1njDQZRb5d9F86FNUKcEjQQLUEMx3oFVKCYU7GvOwjR0P/uh0eYF9Lx2Rc9rUzTy2OsXtCHr+dKyZJdQks0sZPGvFXcTBhI1Ig0qjdtRh4tyMVLfh6LBaDjSzm5DNiSbNmTHcr33ZwgjXlvQe5mPPBhjUYTmlcgVaTTdXZBQnxP83Y7a1R1ICEUHLOPpypuYVi1rEBpHs7ABuu8idisjGJ4rcGwst1Pewc9EHo1omPAFlIdrjkNS3GYi9csz2ItwtBuIojuL8A5awPPHuFfaLaHO+M2C9deCMykeHlsPHkCGvyKUguNr5VEBDTyunJ+7qVzfYCPU2ap9n4ORPWoKElhJiU5H406pTmIyya71JI5C91iWo1FN2vUPkTFnLeQHBo2ei5GPehMUfLQ8MgZNQR6NIwi1gihXI3VyQBXXFMcSQmma17UZZT4KzDkGGSBtWlEsx5VVnKMavksYnFWUMdR2hDwBtyaQha3QSF+tfSnKtWhaYVyLXdBgtRrSGqKYQLG1Cp5rdTRQTEb9Yg6yq9yPu52PRVNqO6U6jbPROzu6wDUORrODNOPqUmUQMkDVOpx4PGEEYAv5U4ANKyFhZxvUQHPmG1h6Lqu90LSvSHaeoTf1+ebAJMKlzrKyKnqermddLdsj4T0lOMd04hdRfZFiXqE90UDp0nZuptKA24ymLlk+ZONiKDKG581LuQ933NIXjlaKS+g41iDssA3UJ318K5be579GIZ980UYFmqc+Sv2WdtstZ/mdcK978XkzE9mL8gxW2yEPiCvhD2S3akdBUYbDkZCqhlvJ7gYH9YVPqLQDeToB6yLD0boFji2hRnw9YVam/3xZOd9Hz2WdjOW7IcNhUrbv+KDOaDbxrcR/e7LZOr8JWLM5nnwfxB0WXEfFQOkXg+3YrIgs7g0o38Ae1ZIyDxuRlhCdeswmOf5iWeQmFNw3ktDmkMQBSDtNShU37k8TTVlCbnmXq39zNKLPR/aHR9C0dX1kLI26ce8n+cO2Njsij+Bce4cXDB2bI5F7zQgAexn5uDDwdmQMs8Nlq8nr76wsQh1uFNks//uhzjYjocymwa9JCmshXNg3ygBkUDb5HU8gV+sZyIC/JeWCYRGKo+lJNmPiRJTV6lpWwOPxpDCR9EQukCB+B7k+42wSJTT1iwbIDUbCeiWrbCthgFgzsgfsjYzEB2ElPiFNop1s7uH1iJlGeDyebDShzpw0PTB8iuwFcZgP2kZDslfGLRiimC9nJ7mSBwRlsqwFMo3KSE+P5/+sgzJgowxzFVzGMd+QTHMH3o7Ue6Mx7EOYY9ID5d7YcRldUYe21wptIFwC4DTClbFANg97ab2hQT1pmk0flOC2bUo5zzLKHigK9D3CaLsN0KhXNFqxs9INhXLba3zamHU3N0aBXRNRYN8WKH7hUdwxE89SqdafHNTVAxk+TeBcV+T1sAMUxyL7Rppr9QJStAX7yzeezsUWKGt1dPCvH+WW9YXIqzEaNcA5aA57K8VTvTsrn6LOeC7qoHGLtzwT7DuOMPKxBa2k9Uf0gSCXYfDLKCcpGhW8NXovrWjJwC2C+iegBZHsVbn2QULjsoT7GBbcwzjis4M9HZzeKHltBDImHUjoadot2Jcl7P1S1IBPxquXaZyIhGkag5FQnoLSAH4b2XeYo/zBVH57tBEJbRO41gtl9sYN6DNJTifvgjwhRZfz93QAmtHiOH3Q6H4Y8nGbNRlPyVhP3+D4hyj2UZhljS4or+e4HMesh7SIaSjBzdV5u6Nl8voXvK7+SANIWl7ul5QLqFiWxhLnntowABnCXkQx9iOQ6vogGll+g1aiirIQNWrDWKR+fooiHnes7yV3GhrRKtRnk82y34C8GoNQjsS2uNPIj0Rh4FkWkbE5Hdkp4hYO3hvFT0wh2xqWng7OaoTf1BhIuFR+Us5DEwqfNVFy30BzWE/96IamBcullLuY/Gnm2+PO9vUsg2yHjEw/IFQPpxJqgYdRuRy64eeUR+aNRsEzu5FudffUlxL53cXDqbH2770SHZvhKMruM5SYcxWhpXke+tzZFsi/3R9NMV5FdokehCnWL6F1JJpQEk7cil2ez4dX04tUVd7j8Xg8Ho/H4/F4PB6P53Phf+i3i7/M0C0DAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "716d4d86",
   "metadata": {},
   "source": [
    "# Kullback-Leibler Divergence (KL)\n",
    "\n",
    "KL divergence score, describes how much one probability distribution differs from a different probability distribution.\n",
    "\n",
    "It is fairly widely used in the data mining literature and as a key ingredient in Variational Auto Encoders. The concept was originated in probability theory and information theory.\n",
    "\n",
    "A KL of zero reflects two distributions that are the same\n",
    "\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "- Which is easier to read and maintain?\n",
    "\n",
    "- Which is faster?\n",
    "\n",
    "### Naive Loop\n",
    "\n",
    "```python\n",
    "\n",
    "def kl_divergence(p, q):\n",
    "    s = []\n",
    "    for i in range(len(p)):\n",
    "        s.append(p[i] * log(p[i]/q[i]))\n",
    "    return sum(s)\n",
    "```\n",
    "\n",
    "### NumPy vector\n",
    "\n",
    "```python\n",
    "def kl_divergence(p, q):\n",
    "    return np.sum(p * np.log(p/q))\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a0460e3",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def kl_divergence_slow(p, q):\n",
    "    s = []\n",
    "    for i in range(len(p)):\n",
    "        s.append(p[i] * log(p[i]/q[i]))\n",
    "    return sum(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "279bb501",
   "metadata": {},
   "source": [
    "# Compute, Print Naive KL Loop, Naive Softmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e96c78d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import log\n",
    "timing = {}\n",
    "start = time.time()\n",
    "np.random.seed(seed=42)\n",
    "BIG = 10_000_000\n",
    "p = softmax_slow(np.random.rand(BIG))\n",
    "q = softmax_slow(np.random.rand(BIG))\n",
    "# calculate (P || Q)\n",
    "kl_pq = kl_divergence_slow(p, q)\n",
    "print('KL(P || Q): %.3f nats' % kl_pq)\n",
    "# calculate (Q || P)\n",
    "kl_qp = kl_divergence_slow(q, p)\n",
    "print('KL(Q || P): %.3f nats' % kl_qp)\n",
    "timing['naiveLoop'] = time.time() - start\n",
    "\n",
    "# plot of distributions\n",
    "from matplotlib import pyplot\n",
    "\n",
    "print(f'timing Naive KL w Softmax Loops: {timing[\"naiveLoop\"]:.2f} sec')\n",
    "print('P=%.3f Q=%.3f' % (sum(p), sum(q)))\n",
    "del p\n",
    "del q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4805e510",
   "metadata": {},
   "source": [
    "# SciPy equivalent: rel_entr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e10f2d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# example of calculating the kl divergence (relative entropy) with scipy\n",
    "from scipy.special import rel_entr\n",
    "# define distributions\n",
    "start = time.time()\n",
    "np.random.seed(seed=42)\n",
    "p = softmax_numpy(np.random.rand(BIG))\n",
    "q = softmax_numpy(np.random.rand(BIG))\n",
    "# calculate (P || Q)\n",
    "\n",
    "###### SciPy equvalent ##################\n",
    "kl_pq = rel_entr(p, q)\n",
    "#########################################\n",
    "\n",
    "print('KL(P || Q): %.3f nats' % sum(kl_pq))\n",
    "# calculate (Q || P)\n",
    "\n",
    "###### SciPy equvalent ##################\n",
    "kl_qp = rel_entr(q, p)\n",
    "#########################################\n",
    "\n",
    "print('KL(Q || P): %.3f nats' % sum(kl_qp))\n",
    "timing['SciPy'] = time.time() - start\n",
    "print(f'timing SciPy: {timing[\"SciPy\"]:.2f} sec')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e100dace",
   "metadata": {},
   "source": [
    "# Compute Naive KL Loop, NumPy Softmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c081108d",
   "metadata": {},
   "outputs": [],
   "source": [
    "start = time.time()\n",
    "np.random.seed(seed=42)\n",
    "p = softmax_numpy(np.random.rand(BIG))\n",
    "q = softmax_numpy(np.random.rand(BIG))\n",
    "# calculate (P || Q)\n",
    "kl_pq = kl_divergence_slow(p, q)\n",
    "print('KL(P || Q): %.3f nats' % kl_pq)\n",
    "# calculate (Q || P)\n",
    "kl_qp = kl_divergence_slow(q, p)\n",
    "print('KL(Q || P): %.3f nats' % kl_qp)\n",
    "timing['naiveLoopFastSoftMax'] = time.time() - start\n",
    "\n",
    "# plot of distributions\n",
    "from matplotlib import pyplot\n",
    "# define distributions\n",
    "print(f'timing NaiveLoopNumPySoftmax: {timing[\"naiveLoopFastSoftMax\"]:.2f} sec')\n",
    "print('P=%.3f Q=%.3f' % (sum(p), sum(q)))\n",
    "del p,q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8ec8ad3",
   "metadata": {},
   "source": [
    "# Compute NumPy KL Loop, NumPy Softmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6d817e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def kl_divergence(p, q):\n",
    "    return np.sum(p * np.log(p/q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e86c7e4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define distributions\n",
    "start = time.time()\n",
    "np.random.seed(seed=42)\n",
    "p = softmax_numpy(np.random.rand(BIG))\n",
    "q = softmax_numpy(np.random.rand(BIG))\n",
    "# calculate (P || Q)\n",
    "kl_pq = kl_divergence(p, q)\n",
    "print('KL(P || Q): %.3f nats' % np.sum(kl_pq))\n",
    "# calculate (Q || P)\n",
    "kl_qp = kl_divergence(q, p)\n",
    "print('KL(Q || P): %.3f nats' % np.sum(kl_qp))\n",
    "timing['NumPyFastAll'] = time.time() - start\n",
    "print(f'timing naiveFastAll: {timing[\"NumPyFastAll\"]:.2f} sec')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c004863-8a6b-4972-857a-1e815c52c75f",
   "metadata": {},
   "source": [
    "# Compute NumPy KL Loop, PyTorch Softmax\n",
    "\n",
    "Something is wrong with PyTorch code \n",
    "\n",
    "See if you can get it to match NumPy code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b38659e3-fe22-4655-9863-568c63c9312b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch.nn.functional as F\n",
    "start = time.time()\n",
    "\n",
    "kl_loss = torch.nn.KLDivLoss(reduction=\"batchmean\")\n",
    "# input should be a distribution in the log space\n",
    "input = torch.from_numpy(p)\n",
    "target = torch.from_numpy(q)\n",
    "# Sample a batch of distributions. Usually this would come from the dataset\n",
    "output = kl_loss(input, target)\n",
    "print(output)\n",
    "timing['pytorch_NOT_SAME_RESULT'] = time.time() - start\n",
    "\n",
    "kl_loss = torch.nn.KLDivLoss(reduction=\"batchmean\", log_target=True)\n",
    "log_target = F.log_softmax(input, dim=0)\n",
    "output = kl_loss(input, log_target)\n",
    "print(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9c5d5d5-615f-456e-90c2-582a215aba44",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,6))\n",
    "plt.title(\"Time taken to process\" ,fontsize=12)\n",
    "plt.ylabel(\"Time in seconds\",fontsize=12)\n",
    "plt.yscale('log')\n",
    "plt.xlabel(\"Various types of operations\",fontsize=14)\n",
    "plt.grid(True)\n",
    "plt.xticks(rotation=-60)\n",
    "plt.bar(x = list(timing.keys()), height= list(timing.values()), align='center',tick_label=list(timing.keys()))\n",
    "print('Acceleration : {:5.0f} X'.format(timing['naiveLoop']/(timing['pytorch_NOT_SAME_RESULT'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fd6081e-8569-4389-9cca-78df26610067",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "print (torch.__version__)\n",
    "\n",
    "m = torch.nn.Softmax (dim=1)\n",
    "in_tensor_before_softmax = torch.Tensor ([0.1, 0.2, 0.4, 0.3])\n",
    "in_tensor = m (in_tensor_before_softmax.view (-1,4))\n",
    "out_tensor_before_softmax = torch.Tensor ([0.7, 0.1, 0.1, 0.1])\n",
    "out_tensor = m (out_tensor_before_softmax.view (-1,4))\n",
    "\n",
    "import torch.nn.functional as F\n",
    "kl_loss = torch.nn.KLDivLoss (reduction = 'batchmean')\n",
    "loss = torch.nn.CrossEntropyLoss()\n",
    "kl_output = kl_loss (input = F.log_softmax (in_tensor_before_softmax, dim=-1), target = out_tensor)\n",
    "cross_ent = loss (input = in_tensor_before_softmax.view (-1,4), target = out_tensor)\n",
    "ent = loss (input=out_tensor_before_softmax.view (-1,4), target = out_tensor)\n",
    "kl_output_using_ce = cross_ent - ent\n",
    "print (kl_output, kl_output_using_ce)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0164a13f-035b-4ed3-8c0e-68fd65e070c9",
   "metadata": {},
   "source": [
    "# Experiment to get same KM_output as NumPy\n",
    "\n",
    "Experiment not working yet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ce993ff-96bf-4d94-95bb-e5e4f861d869",
   "metadata": {},
   "outputs": [],
   "source": [
    "kl_loss = torch.nn.KLDivLoss (reduction = 'none')\n",
    "loss = torch.nn.CrossEntropyLoss (reduction = 'none')\n",
    "kl_output_unr = kl_loss (input = F.log_softmax (in_tensor_before_softmax, dim=-1), target = out_tensor)\n",
    "cross_ent_unr = loss (input = in_tensor_before_softmax.view (-1,4), target = out_tensor)\n",
    "ent_unr = loss (input=out_tensor_before_softmax.view (-1,4), target = out_tensor)\n",
    "kl_output_using_ce_unr = cross_ent_unr - ent_unr\n",
    "print (kl_output_unr)\n",
    "\n",
    "print (kl_output_using_ce_unr)\n",
    "\n",
    "print (kl_output_unr.mean(), kl_output_unr.sum(), kl_output_using_ce_unr)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "488819a9-3830-49f2-8f6e-c8a09196fecc",
   "metadata": {},
   "source": [
    "# Notices and Disclaimers\n",
    "\n",
    "Intel technologies may require enabled hardware, software or service activation.\n",
    "No product or component can be absolutely secure. \n",
    "\n",
    "Your costs and results may vary. \n",
    "\n",
    "© Intel Corporation. Intel, the Intel logo, and other Intel marks are trademarks of Intel Corporation or its subsidiaries. Other names and brands may be claimed as the property of others. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a30d52ef",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PyTorch (AI kit)",
   "language": "python",
   "name": "c009-intel_distribution_of_python_3_oneapi-beta05-pytorch"
  },
  "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}