gaussian_evaluation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Investigating speed of scipy.stats.norm vs custom gaussian function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "\n",
    "def gaussian(xvals, sigma, centre=0):\n",
    "    # evaluate gaussian function with defined sigma and center at x\n",
    "    gauss_eval = np.exp(-0.5 * ((xvals - centre) / sigma)**2) \\\n",
    "                    / (sigma * np.sqrt(2 * np.pi))\n",
    "    return gauss_eval\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, we can time both functions with arguments similar in size to those used within `fast_broaden.find_coeffs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "xvals = np.linspace(0,100,101)\n",
    "sigma = np.linspace(1, 5, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17.2 µs ± 436 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit gaussian(xvals, sigma[:,np.newaxis])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "81.7 µs ± 565 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit norm.pdf(xvals, scale=sigma[:,np.newaxis])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For these inputs, the custom gaussian function is the faster option. Next, we can test how the speed of the functions scales when the length of the xvals array increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAU/klEQVR4nO3df6xf9X3f8edr1zFJSAPEtSqC0XCFW2SmtaFXlLtW0VW8BpNWuH8g1WRdUceG1MKWLqsqW1WmjfyBmCbSViVprUDH0qSGudlqRWpoY+I/Jt0ZrpukwXa93AEdZmTcEEK2aMO1894f32N2uVzu/dq+3x/ne54P6avv+X7O53vu56P71Xmd8zm/UlVIkrrnb426AZKk0TAAJKmjDABJ6igDQJI6ygCQpI4yACSpo/oKgCQ7k5xMspBkzwrzL0nyaDP/SJJrmvJNSb6c5H8n+d1l3/mJJF9vvvM7SbIuPZIk9WXNAEgyBTwI3AJsB25Psn1ZtTuBV6rqWuATwP1N+f8FPgb8+gqL/hTwT4BtzWvnhXRAknRh+tkDuBFYqKpnquo0sB/YtazOLuCRZvoAsCNJqup7VfWf6QXB65JcCby7qv5L9a5E+/fAz19EPyRJ52lDH3WuAp5f8vkU8JNvVaeqziR5FdgEfGuVZZ5atsyrVqqY5C7gLoBLL730J6677ro+mixJAjh69Oi3qmrzSvP6CYCRqqp9wD6A6enpmp+fH3GLJKk9kvz1W83rZwjoBeDqJZ+3NGUr1kmyAbgMeHmNZW5ZY5mSpAHqJwCeArYl2ZpkI7AbOLiszkHgjmb6NuCJWuUuc1X1IvDdJDc1Z//8EvAn5916SdIFW3MIqBnTvwd4HJgCHq6qY0nuBear6iDwEPCZJAvAt+mFBABJngPeDWxM8vPAB6vqOPCrwL8D3gH8afOSJA1J2nQ7aI8BSNL5SXK0qqZXmueVwJLUUQaAJHWUASBJY2xuDu67r/e+3sb+OgBJ6qq5OdixA06fho0b4dAhmJlZv+W7ByBJY+rw4d7K/+zZ3vvhw+u7fANAksbU7Gxvy39qqvc+O7u+y3cISJLG1MxMb9jn8OHeyn89h3/AAJCksTYzs/4r/nMcApKkjjIAJKmjDABJ6igDQJI6ygCQpI4yACSpowwASeooA0CSLtAgb9Q2DF4IJkkXYNA3ahsG9wAk6QIM+kZtw2AASNIFGPSN2obBISBJugCDvlHbMBgAknSBBnmjtmFwCEiSOsoAkKSOMgAkqaMMAEnqKANAkjrKAJCkjjIAJKmjDABJ6igDQNLEaftdOofFK4ElTZRJuEvnsLgHIGmiTMJdOofFAJA0USbhLp3D4hCQpIkyCXfpHJa+9gCS7ExyMslCkj0rzL8kyaPN/CNJrlkyb29TfjLJzUvK/3mSY0meTvJHSd6+Lj2S1HkzM7B3ryv/tawZAEmmgAeBW4DtwO1Jti+rdifwSlVdC3wCuL/57nZgN3A9sBP4ZJKpJFcB/wyYrqq/A0w19SRJQ9LPHsCNwEJVPVNVp4H9wK5ldXYBjzTTB4AdSdKU76+q16rqWWChWR70hp/ekWQD8E7gf1xcVyRJ56OfALgKeH7J51NN2Yp1quoM8Cqw6a2+W1UvAP8W+O/Ai8CrVfVnK/3xJHclmU8yv7i42EdzJUn9GMlZQEmuoLd3sBV4L3Bpkl9cqW5V7auq6aqa3rx58zCbKUkTrZ8AeAG4esnnLU3ZinWaIZ3LgJdX+e7fB56tqsWq+hvg88Dfu5AOSGoXr9IdH/2cBvoUsC3JVnor793Ah5fVOQjcAcwBtwFPVFUlOQh8LskD9Lb0twFPAt8HbkryTuD/ADuA+XXoj6Qx5lW642XNPYBmTP8e4HHgBPBYVR1Lcm+SW5tqDwGbkiwAHwX2NN89BjwGHAe+CNxdVWer6gi9g8V/AXy9ace+de2ZpLHjVbrjJVU16jb0bXp6uubn3VGQ2so9gOFLcrSqplea55XAkobGq3THiwEgaahmZlzxjwtvBidJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASDpdXNzcN99vXdNvg2jboCk8TA3Bzt2wOnTsHEjHDoEMzOjbpUGyT0ASQAcPtxb+Z8923s/fHjULdKgGQCSAJid7W35T0313mdnR90iDZpDQJKA3nDPoUO9Lf/ZWYd/usAAkPS6mRlX/F3iEJAkdVRfAZBkZ5KTSRaS7Flh/iVJHm3mH0lyzZJ5e5vyk0luXlJ+eZIDSf4qyYkkbndI0hCtGQBJpoAHgVuA7cDtSbYvq3Yn8EpVXQt8Ari/+e52YDdwPbAT+GSzPIDfBr5YVdcBPwacuPjuSJL61c8ewI3AQlU9U1Wngf3ArmV1dgGPNNMHgB1J0pTvr6rXqupZYAG4McllwPuBhwCq6nRVfeeieyNJ6ls/AXAV8PySz6eashXrVNUZ4FVg0yrf3QosAn+Q5CtJPp3k0pX+eJK7kswnmV9cXOyjuZKkfozqIPAG4AbgU1X1PuB7wJuOLQBU1b6qmq6q6c2bNw+zjZI00foJgBeAq5d83tKUrVgnyQbgMuDlVb57CjhVVUea8gP0AkGSNCT9BMBTwLYkW5NspHdQ9+CyOgeBO5rp24Anqqqa8t3NWUJbgW3Ak1X1TeD5JD/afGcHcPwi+yJJOg9rXghWVWeS3AM8DkwBD1fVsST3AvNVdZDewdzPJFkAvk0vJGjqPUZv5X4GuLuqzjaL/qfAZ5tQeQb45XXumyRpFeltqLfD9PR0zc/Pj7oZktQaSY5W1fRK87wSWGoB79OvQfBeQNKY8z79GhT3AKQx5336NSgGgDTmvE+/BsUhIGnMeZ9+DYoBILWA9+nXIDgEJEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoB0kXxYi9rKm8FJF8GHtajN3AOQLoIPa1GbGQDSRfBhLWozh4Cki+DDWtRmBoB0kXxYi9rKISBJ6igDQJLG2QDPM3YISJLG1YDPM3YPQJLG1YDPMzYAJGlcDfg8Y4eAJGlcDfg8YwNAksbZAM8zdghIkjrKAJCkjjIAJKmjDABJulAtfxhEXwGQZGeSk0kWkuxZYf4lSR5t5h9Jcs2SeXub8pNJbl72vakkX0nyhYvuiSSdM4wV87mLtD72sd57C0NgzbOAkkwBDwI/A5wCnkpysKqOL6l2J/BKVV2bZDdwP/ALSbYDu4HrgfcCX0ryI1V1tvneR4ATwLvXrUfSEnNz3qmzc4b1lJ6VLtJq2Y+snz2AG4GFqnqmqk4D+4Fdy+rsAh5ppg8AO5KkKd9fVa9V1bPAQrM8kmwBfhb49MV3Q3qzCdhA04UY1lN6JuBhEP0EwFXA80s+n2rKVqxTVWeAV4FNa3z3t4DfAL6/2h9PcleS+STzi4uLfTRX6vFpXR01rBXzuYu0Pv7x1j4LdCQXgiX5OeClqjqaZHa1ulW1D9gHMD09XYNvnSbFufXAuZGAFm6g6UIM8yk9LX8YRD8B8AJw9ZLPW5qyleqcSrIBuAx4eZXv3grcmuRDwNuBdyf5w6r6xQvqhbQCn9bVYS1fMQ9LPwHwFLAtyVZ6K+/dwIeX1TkI3AHMAbcBT1RVJTkIfC7JA/QOAm8DnqyqOWAvQLMH8Ouu/DUIrgekt7ZmAFTVmST3AI8DU8DDVXUsyb3AfFUdBB4CPpNkAfg2vZCgqfcYcBw4A9y95AwgSdIIpao9w+rT09M1Pz8/6mZIUmskOVpV0yvN80pgSeooA0CSOsoA0Ei0/BYq0kTwgTAaumFdqS9pde4BaOi8QlcaDwaAhm4CbqEiTQSHgDR0XqErjQcDQCPhFbrS6DkEJEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHWUASBJHWUASFJHGQCS1FEGgCR1lAEgSR1lAEhSRxkAktRRBoAkdZQBIEkd1VcAJNmZ5GSShSR7Vph/SZJHm/lHklyzZN7epvxkkpubsquTfDnJ8STHknxk3XokSerLmgGQZAp4ELgF2A7cnmT7smp3Aq9U1bXAJ4D7m+9uB3YD1wM7gU82yzsD/Iuq2g7cBNy9wjIlTaK5Objvvt67RmpDH3VuBBaq6hmAJPuBXcDxJXV2Af+qmT4A/G6SNOX7q+o14NkkC8CNVTUHvAhQVf8ryQngqmXLlDRp5uZgxw44fRo2boRDh2BmZtSt6qx+hoCuAp5f8vlUU7Zinao6A7wKbOrnu81w0fuAIyv98SR3JZlPMr+4uNhHcyWNrcOHeyv/s2d774cPj7pFnTbSg8BJ3gX8MfBrVfXdlepU1b6qmq6q6c2bNw+3gZLW1+xsb8t/aqr3Pjs76hZ1Wj9DQC8AVy/5vKUpW6nOqSQbgMuAl1f7bpK30Vv5f7aqPn9BrZe0vubmelvls7ODGZqZmekN+wzyb6hv/QTAU8C2JFvprbx3Ax9eVucgcAcwB9wGPFFVleQg8LkkDwDvBbYBTzbHBx4CTlTVA+vTFUkXZVjj8zMzrvjHxJpDQM2Y/j3A48AJ4LGqOpbk3iS3NtUeAjY1B3k/CuxpvnsMeIzewd0vAndX1Vngp4B/CHwgyVeb14fWuW+Szofj852Tqhp1G/o2PT1d8/Pzo26GNJk8Q2ciJTlaVdMrzetnCEjSqA16bB4cn+8gA0Aad8PcMnd8vlO8F5A07hyb14AYAHoTr9QfM547rwFxCEhv4HHAMeTYvAbEANAbrDTa4PpmDDg2rwFwCEhv4GiD1B3uAegNHG2QusMA0JtM1GjDMM6fl1rKAGgZ12fnYVhHtP2nqKUMgBbxDJ3zNIwj2v5T1GIeBG4Rrwc6T8M4ou0/RS3mHkCLnFufndvY9AydNQzjiLb/FLWYdwNtGYebx5D/FI2x1e4GagBI0gRbLQA8BqDR8IZD0sh5DEDD55kz0lhwD2CduEF7HjxzRhoL7gGsAzdoz5NnzkhjwT2AdTBxG7SD3p05d3rmxz9uWkoj5B7AOpioDdph7c5M1A2HpHZyD2AdDHWDdtBb5xO3OyPprbgHsE6GskE7jK3zidqdkbQaA6BNhnFzMx8IIHWGAdAmw9o6d3xe6oROBMDE3KrFrXNJ62jiA2DiztF361zSOpn4s4CGdlKLlwJLapmJ3wMYyrD5xO1mSOqCiQ+AoQybD+PsHElaZxMfAAAzzDHDYWAW8KlQkgRdCIBhDM94do6kFpr8ABjW8Ixn50hqmYk/C+j14ZmpKYdnJGmJvgIgyc4kJ5MsJNmzwvxLkjzazD+S5Jol8/Y25SeT3NzvMteNtx6WpBWtOQSUZAp4EPgZ4BTwVJKDVXV8SbU7gVeq6toku4H7gV9Ish3YDVwPvBf4UpIfab6z1jLXj8MzkvQm/ewB3AgsVNUzVXUa2A/sWlZnF/BIM30A2JEkTfn+qnqtqp4FFprl9bNMSdIA9RMAVwHPL/l8qilbsU5VnQFeBTat8t1+lglAkruSzCeZX1xc7KO5kqR+jP1B4KraV1XTVTW9efPmUTdHkiZGPwHwAnD1ks9bmrIV6yTZAFwGvLzKd/tZpiRpgPoJgKeAbUm2JtlI76DuwWV1DgJ3NNO3AU9UVTXlu5uzhLYC24An+1ymJGmA1jwLqKrOJLkHeByYAh6uqmNJ7gXmq+og8BDwmSQLwLfprdBp6j0GHAfOAHdX1VmAlZa5/t2TJL2V9DbU2yHJIvDXo27HKn4Q+NaoG7FO7Mv4mZR+gH0Zpr9dVSseQG1VAIy7JPNVNT3qdqwH+zJ+JqUfYF/GxdifBSRJGgwDQJI6ygBYX/tG3YB1ZF/Gz6T0A+zLWPAYgCR1lHsAktRRBoAkdZQBsIYkDyd5KcnTS8rek+TPk3yjeb+iKU+S32mecfCXSW5Y8p07mvrfSHLHSn9rwP24OsmXkxxPcizJR1rcl7cneTLJ15q+/OumfGvzPIqF5vkUG5vy835exZD7M5XkK0m+0PJ+PJfk60m+mmS+KWvd76tpw+VJDiT5qyQnksy0tS+rqipfq7yA9wM3AE8vKfs3wJ5meg9wfzP9IeBPgQA3AUea8vcAzzTvVzTTVwy5H1cCNzTTPwD8V2B7S/sS4F3N9NuAI00bHwN2N+W/B/xKM/2rwO8107uBR5vp7cDXgEuArcB/A6ZG8Bv7KPA54AvN57b24zngB5eVte731bTjEeAfN9Mbgcvb2pdV+znqBrThBVzDGwPgJHBlM30lcLKZ/n3g9uX1gNuB319S/oZ6I+rTn9B7IE+r+wK8E/gL4CfpXY25oSmfAR5vph8HZprpDU29AHuBvUuW9Xq9IbZ/C3AI+ADwhaZdretH83ef480B0LrfF72bWT5Lc5JMm/uy1sshoAvzQ1X1YjP9TeCHmumLfv7BMDRDB++jt+Xcyr40wyZfBV4C/pzeVu93qvc8iuXtOt/nVQzTbwG/AXy/+byJdvYDoIA/S3I0yV1NWRt/X1uBReAPmqG5Tye5lHb2ZVUGwEWqXrS35lzaJO8C/hj4tar67tJ5bepLVZ2tqh+ntwV9I3DdaFt0/pL8HPBSVR0ddVvWyU9X1Q3ALcDdSd6/dGaLfl8b6A37fqqq3gd8j96Qz+ta1JdVGQAX5n8muRKgeX+pKR/r5x8keRu9lf9nq+rzTXEr+3JOVX0H+DK9oZLL03sexfJ2ne/zKoblp4BbkzxH77GoHwB+m/b1A4CqeqF5fwn4j/SCuY2/r1PAqao60nw+QC8Q2tiXVRkAF2bp8w/uoDeefq78l5qzAm4CXm12GR8HPpjkiubMgQ82ZUOTJPRu232iqh5YMquNfdmc5PJm+h30jmWcoBcEtzXVlvflfJ5XMRRVtbeqtlTVNfQO6j5RVf+AlvUDIMmlSX7g3DS938XTtPD3VVXfBJ5P8qNN0Q56t7RvXV/WNOqDEOP+Av4IeBH4G3pbBnfSG3c9BHwD+BLwnqZugAfpjUd/HZhespx/BCw0r18eQT9+mt4u618CX21eH2ppX/4u8JWmL08D/7Ip/2F6K74F4D8AlzTlb28+LzTzf3jJsn6z6eNJ4JYR/s5m+f9nAbWuH02bv9a8jgG/2ZS37vfVtOHHgfnmN/af6J3F08q+rPbyVhCS1FEOAUlSRxkAktRRBoAkdZQBIEkdZQBIUkcZAJLUUQaAJHXU/wPiWyvdTFweNwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import timeit\n",
    "\n",
    "sigma = np.linspace(1, 10, 30)\n",
    "lengths = np.arange(500,7000,500)\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "for l in lengths:\n",
    "    xvals = np.linspace(0, 1000, l)\n",
    "\n",
    "    def run_gaussian():\n",
    "        gaussian(xvals, sigma[:,np.newaxis])\n",
    "\n",
    "    def run_norm_pdf():\n",
    "        norm.pdf(xvals, scale=sigma[:,np.newaxis])\n",
    "\n",
    "    run_time_gauss = timeit.Timer(stmt=run_gaussian).repeat(5,1)\n",
    "    run_time_norm = timeit.Timer(stmt=run_norm_pdf).repeat(5,1)\n",
    "\n",
    "    plt.plot(l, np.mean(run_time_gauss),'r.')\n",
    "    plt.plot(l, np.mean(run_time_norm),'b.')    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regardless of the length of the xvals array, the custom gaussian function is faster. Finally, it should be checked that the two functions produce the same output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "xvals = np.linspace(-5,5,101)\n",
    "sigma = 2\n",
    "scipy_norm = norm.pdf(xvals, scale=sigma)\n",
    "gaussian_eval = gaussian(xvals, sigma)\n",
    "\n",
    "test = gaussian_eval.all() == scipy_norm.all()\n",
    "\n",
    "print(test)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}