Speed comparison between models in Sherpa, Astropy and astromodels

sherpa_astropy_astromodels.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:sherpa.astro.all:failed to import sherpa.astro.xspec; XSPEC models will not be available\n"
     ]
    }
   ],
   "source": [
    "from sherpa.astro import ui as sherpa\n",
    "\n",
    "import logging\n",
    "logging.getLogger('sherpa').propagate = 0\n",
    "\n",
    "import astromodels as am\n",
    "\n",
    "import astropy.modeling.models as apm\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First define a power-law implemented with numpy. This is the reference, as anything you add around this will slow things down:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reference implementation, i.e., a minimal implementation\n",
    "# of a power law. Whatever we do, we will never be faster\n",
    "# than this if we use pure python (no cython)\n",
    "\n",
    "def plain_numpy_powerlaw(energies, norm, index):\n",
    "    \n",
    "    return norm * np.power(energies, index)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define the same power law in all 3 packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sherpa_powerlaw = sherpa.powerlaw(\"test\")\n",
    "sherpa_powerlaw.index = -2\n",
    "sherpa_powerlaw.refer = 1.0\n",
    "sherpa_powerlaw.ampl = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "astromodels_powerlaw = am.powerlaw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "astropy_powerlaw = apm.PowerLaw1D()\n",
    "\n",
    "# Astropy uses the formula with the positive index\n",
    "\n",
    "astropy_powerlaw.alpha = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluate the powerlaws on an array of energies and plot them, just to check that they are indeed identical:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGBJJREFUeJzt3X2QXWV9wPFvIq8FWRCwrRSziCRxFBsEFKWwl91Naoep\ntmgYSRqFUpXSTaLIRKcWQ9pOZ6pACqEVadxSKMSCBkekLTYiyJCgEMKQkkSJCqVFqW00rNagDds/\nzlmy92Zf7st57nNevp+ZO7v33nPP87szm/PL7zxvIEmSJEmSJEmSJEmSJEmSJEmSJEmScudlsQOY\nwsHAp4DfAfYCO+OGI0nKu37ggvT3G2MGIknaZ2aX2xsGngO2Nrx+NrAdeBJYmr52MvCd9PdDuxKd\nJCl3zgJOYf/EsYUkecwCdgDHAOcA70nf/0y3ApQkTa3bFccDwI8aXutJf34deBr4CvAWYCNwGrAG\nWN+tACVJUzsgdgDA6SRVxphtwBnA3cDlU33wkMNf/vM9Pxk5MGBsklRG3wFe2+6Hu11xZGrPT0YO\n7D9/6e+dtvqqF+cMr2Vw0dA1wIyCPlaVpM1Oz9nu51v5XLPHTndcp+8X5RHre4Rot9NztvP5Vj7T\n7LHNHDfVMSdSML3U93H0kPRxjFkDnNvkuUaBK5kxo7bg0uXbj75z/ejAihW7+s9bflQ2oXZVrSRt\ndnrOdj/fyueaPXa646Z7/8om28m7Wona7fSc7Xy+lc80e2wzx010TI3k73K0yXZyo5fJO8d72dc5\n3oy6L99//tLFp62+au+c4bWjg4uGru4wTqlTV8YOQJpEoRLHOuBZ4AXgGeCi9PU+kuG4O4FlLZxv\nvy/ft3zljAWXLt9R8OpD5VCLHYA0iUIljqwlt6om+Adq9SFJ+6mRwa2qGVlEEtEoU3yHvuUrZxz8\nix9v3zy/b868TQ/9aHTnCyfeu/7axuHAklQ1U147p1PqxDGm//yli59/66ybR3qOnHn8hseu2XDb\n9R/pQmySlFeVTxyrgPvSx6Qaq48Xd+6Z87X11/2wCzFKUl7U0sdKKp44WvoOY9XH8z1Hznz1Vx/7\nqw23Xv/hQLFJUl5VvuJo+Tuk1ce2zfP75lp9SKogE0e7H7b6kFRRlU8cTfVxTMa+D0kVUsM+js6y\n5nhWH5IqpPIVR2bfwepDUkWYOLI+qdWHpJIzcYQ4sdWHpBILdu0sgknXqsrK2JpXs4fXjg4uHlod\nqh1J6oIarlXVnaxp9SGpZLxV1a3G7PuQVBImjm42aPUhqQRMHDEatvqQVGCVTxwdzRzvhGteSSqY\nGs4cz8eQMqsPSQVT+YojF9+hofrY/eLOPSdZfUjKKRNH7CDGs/qQVAAmjthBNLL6kJRzJo7YQUym\nrvpwr3NJ+WHiiB3EVBrnfYzufOHEe9df+6PYcUmqtNxfO0MKvlZVVsbWvJozvHZ0cNHQ1bHjkVRJ\nNVyrqlhZ01nnknLCW1Wxg2iVI68kRWbiiB1EO6w+JEVk4ogdRCesPiRFYOKIHUSnrD4kdZmJI3YQ\nWbH6kNQlJo7YQWTJ6kNSF5g4YgcRgtWHpIBMHLGDCMXqQ1Igpb52TqcwM8c7MTbrfPbw2tHBxUOr\nY8cjqbBqOHO8OlnT3QYlZchbVbGD6CZX3JWUARNH7CC6zRV3JXXIxBE7iFgGFg4t2X3mCTeNHNEz\n83irD0nNM3HEDiKmtPrYsXl+32yrD0lNMnHEDiIPxvo+RnqOtPqQNB0TR+wg8sK+D0lNMnHEDiJv\n6qoPZ51L2p+JI3YQefRS9bGgNmfexk27X9y55yTnfUhKmThiB5FnrnklaQImjthB5J1rXklqUNpr\n5wnAWuCOKY7paL2VqnHNK0mp0q9VdQewcJL3Sps1Q7H6kEQBrp3DwHPA1obXzwa2A08CS6f4vBVH\nAFYfUqXl/tp5FnAK+yeOLSTJYxawAzgGWAKsBl417jgTRyB9y1fOWHDp8h1H37l+dGDFil395y0/\nKnZMkrqiENfOXuoTRw9J4hhzHXBuw2deAdxAUpF8dJLzFuLL591Y9TFneO3o4KKhq2PHIym4QvRx\n9AJ3ASenzweBi4EL0ueXAMcBV7R43lFg1bjn96UPtchZ51Kp1ajf8G4lOe/jgP0rjkFg3bjnlwB/\n1sZ5rTgyVld92PchlVVH186ZWUXRooeBueOevx54qM1zXUnJt47tpntvX3PrYU+NHDDrka3f2vLu\n/g/Z9yGVSo3kmlkIvUzeOd7Lvs7xVllxBGTfh1Raub92rgOeBV4AngEuSl/vIxmOuxNY1ua5c//l\ni86RV1IpVfraOYq3qrrC6kMqhRrJNbMQo6pCyf3sxzJx1rlUGi5yGDuIqnHFXanwKp84VuH8ja5z\nvw+pkGrpo6N5HGVIHEX/DoVm9SEVUuUrjqJ/h8JLq49tm+f3zZ236SGrDyn/TByxg1DC6kMqjMon\nDvs4csQ1r6Rcq2EfhxVHXvWfv3Tx82/rvXnkiJ6Zx2947JoNt13/kdgxSXpJ5SuOon+H0kqrjx2b\n5/fNtvqQcsXEETsITW1g4dCS3WeecJPVh5QblU8c9nEUgH0fUi7UsI/DiqNoxkZejfQcOfN4R15J\nsQSvOAaBDQ2vvQ/4+3YbzZCJo4Bc80qKLvi18wHg08BhwK+QbAH7hZANtqDSSwMX3diKu7PdbVDq\ntuCr484EPgJ8MG1sJXBbJ41myIqj4BpmnVt9SN0R/Np5NHA7cA/wBPCx0A22wIqjJKw+pK4Kfu38\nNnBx+vsvAWuAjaEbbZIbOZVI426D55y37NjYMUklUyODjZyaMWuC1/pCN9okK44SsvqQggvex9E3\nSSNf76ThjNjHUVLu9yEFFXw47pfZlzheAbyZZLLd/HYbzZCJo+TqVtx11rmUla7PHH8Dyciqhe02\nmiETRwU461zKXNevnQcB27rZ4BTs46iQsb6POcNrRwcXDV0dOx6pwIL3cawZ9/vBwBkkkwA/3knD\nGbHiqBirDykTwW9VXTju958BDwFPt9tgxlzksKLs+5DaUsNFDq04qszqQ2pbsIpj6zSNvrHdRjNk\n4lD9irtWH1IzgiWOucCeKRp4qt1GM2TiEGD1IbUo2LXz0fTnLSFOnhFHValO3cgrZ51Lkwk2quoB\n4LPAJ4DLG44dBdZ30nBGrDi0H/f7kKYV7FbVPOC9JJs2fWmC9y9qt9EMmTg0qbqRV+42KI0XfDju\nHwBr220gMBOHpmT1IU2o60uO5ImJQ02x+pDqmDhiB6FiaNht0BV3VWWVvna6kZNaVrffh2teqVpq\nZLCR01QZ59RpTv7oFO91S6WzptrnvA9VXLBbVfcxdeI4p91GM2TiUEcGFg4t2X3mCTeNHNHjrHNV\niX0csYNQsaXVx47NC2qz523cZPWhKuhK4jgM6AeOGvfaze02miEThzIzsHBoye639d7kmleqgOCJ\n4/3AxcBrgAeBAZL9OBa322iGTBzKlPM+VBHBr50bSXb9eyJ9Phv4SsgGW+BaVQqibuSVa16pfILv\nAPgwcDpwN/AB4AfAvwGv66ThjFhxKBirD5VY8Gvnp0n6Ns4Hvgt8C/iLkA22wIpDwVl9qISCVxzj\nvZwkifx7J41myIpDXWH1oZIJ1jn+OmA78KZJ3ncCoCrHNa9UEsESx9+SjKi6j4nLGicAqpKsPlQC\nwYfjHkKyhex0r4XwTuBc4ADgBuCbDe+bOBSN1YcKLHjieJT9b1dN9FpIrwRWAX/Y8LqJQ1G55pUK\nKti181dJFjrcQZIkTk1/vh3Y1OK5hoHngK0Nr59N0o/yJLB0is9fTbIjYSNHVSkX6vY6d8Vd5V+w\nUVXvAy4ETgMeGff608A/ABtaaOcs4Ccky5ScPO71LcDy9Jz3AL8B/BZJgvoU8H3gL9P3vjrBea04\nlBtWHyqQ4Leq3g18vt0GxuklWapkLHH0kHS8n5I+v44kQdw97jPLSPY9fxh4DPhMwzlHSW5hjbkv\nfUjRjPV9uOaVcqRG/b5FKwmcOA4E3po+Dk4/Mwr8aYtt9VKfOAZJ1sC6IH1+CXAccEUL57TiUC5Z\nfSjnOrp2HtDEMWtILvr3Az9ttyGpSu6/dtUoMLf/v19cvPuts24emXvkrsFDhxx5pcrYBszM4Dy9\n1HeO95D0cYxZQzL0thVuHavc61u+csaCS5fvOPqLd44OrFjx43POW3Zs7JhUWTUy2Dq2GX9NspR6\np3rZf1TVFpKRVb0ko7eOafGcjqpSYbjmlXIk+FpV24C5wH8CPx7X6BtbaGcd0AccDfwX8Ang79LX\nbiDpR7kufbTCPg4VirPOlRPBR1X1TvL6U+02mqGxUVX34WgqFYizzhVJLX0EH1UFyWTAfuBW4Fjg\ncOB77TaaISsOFdZL1ceC2px5GzftfnHnnpOsPtQlwa+dHyD53/y30+fHkWwhmwf2cajw7PtQBMGv\nnQ+QbB07fgTU46EbbZKjqlQKL428unP96MCKFbsceaVAanRpVNWXSeZ7jCWOVwP/FLrRJllxqFSs\nPtQlwUdVLSRZduTNwE3Au0gy1vpOGs6IfRwqHfs+1AVBR1XNAI4nmQD4rvTn54Bn2m0wYyYOlZYj\nrxRQ8MTxOPUr2uaJw3FVaq55pYzVyGA4bjNuJNmJL4/s41AluN+HMha8j2M7MAf4H+AH4xptZeZ4\nKN6qUmVYfShDlZ85buJQpbjfhzIQPHFAvmeO28ehynHNK7WpRpf6OJw5LuWU8z7UpuB9HA+QLKv+\nDfZt8/o49nFIuZBWH9s2z++bO2/TQ877UDOCXzudOS4VgNWHWuDM8dhBSHlh34eaFKxz/BBgT/p7\nL84clwrDWeeaRrDE8SjwJuAWYEm7DQTmqCppEq55pQnUCLyR0wPAZ0m2eb08PXZ03E9vVUkFYPWh\nCQSrOOYB7wXeB3xpgvcvarfRDJk4pCY461wNgl07F6Y/PxDi5BlxVJXUAte8UirYqKrNwKkkw3BP\nmeK4mKw4pBZZfYiAt6puJZklPo+kv6Ox0Xe022iGTBxSm1zzqtKCrlX1RpJO8Isbjh0F7m+30QyZ\nOKQOWH1UVvBr57EhT96hUZLJiLW4YUjFZt9HZdRIrpnB+jiuBZYDd03wnreqpJKx+qiUYLeqTiXp\nIK9N0qi3qqQSsu+jErqyH8fhaUM/bbehQEwcUgANK+665lX5BLt2zgA+BDwL/Ixk3ar/ILl9lZeL\ntfM4pIBccbe0gvVx/D7wfuCTwL+mxy4ALiNZimS4k4YzYsUhBeaaV6UU7Nr5MPCbE7w+CDwSosE2\nWHFIXWL1USrBKo5vA68D9ja8/jJgOzC7k4YzYsUhdZH7fZRGsGvnljbf6yYrDikCq4/CC1Zx7AX+\nd5L3DiXZTjY29+OQIrH6KKQagffjKAJvVUmR1e334byPoujKPI68MnFIOeCs88IxccQOQlLCWeeF\nYeKIHYSkfez7KAQTR+wgJO3Pvc5zzcQROwhJE7P6yC0TR+wgJE3N6iN3TByxg5A0vcYVdx15FZWJ\nI3YQkprnyKtcMHHEDkJSa6w+ojNxxA5CUnvqqg/7PrqptIljLsmmUQcBdwPrJzjGxCEVnCOvoij9\ntfMgYN0k77k6rlQSrrjbVbm/dg4DzwFbG14/m2RfjyeBpZN89h3Ag8DvTvJ+7r+8pOb1LV85Y8Gl\ny7cffef60YEVK3adc96yY2PHVFK5v3aeBZzC/oljC0nymAXsAI4BlgCrgVc1HPulSc6d+y8vqXV1\n1ceioatjx1NCwfbjyFIvcBdwcvq8h2T/jFPS59cB95D0ZYzpA84jifFh4JYJzlv6+3RSVbniblAd\nXTtjbcZ0OkmVMWYbcAb1ieP+9DGdK8f9fh9u6CSVwv3XrhoF5g78cO+S3WeecNPI3J5dg4cMOe+j\nPbX0kYlYFccgcDFwQfr8EuA44IoWz2vFIVVAWn3s2LygNnvexk1WH50rxLWzl/o+jh7q9y1fA5zb\nxnnt45AqZGDh0JLTVl+1d459H50qxLWzl8k7x3vZ1zneqlGSW1W1tiOTVCjpyKsdjrxqS43kmpn7\nxLEOeBZ4AXgGuCh9vY9kOO5OYFmb5879l5cUhvM+OlKIUVWhFOI+naQwXhp5taA2Z97GTbtf3Lnn\nJGedN6XS105vVUmy+mhejQxuVRU941Q6a0rax+qjJaVd5LAZJg5JddxtsCmVTxyrcOKfpHGsPiZV\nSx8rqXjiKPp3kBSI1cekKl9xFP07SArI/T4mZOKIHYSk/LP6qFP5xGEfh6SmWH3YxwFWHJLaYPVh\nxVH07yApgorv92HiiB2EpOIaqz5Geo6cefyGx6qy30flE4d9HJI6UqHqo4Z9HFYckrJTV32Uu++j\n8hVH0b+DpBypyMgrE0fsICSVT8lHXpk4YgchqZxKvOZVpa+d7schKbgS7fdRw/04qp01JXVPyaoP\nb1XFDkJSdZSk78PEETsISdVSgnkfJo7YQUiqpv7zly5+/m29N48c0VO0WecmjthBSKqutPrYsXl+\n3+wCVR8mjthBSNLAwqElu8884aaCVB+VTxyuVSUpFwrQ91HDtaqsOCTlTwHWvKp8xVH07yCphHK+\n5pWJI3YQkjSZnM77MHHEDkKSppLDWecmjthBSFIz6qqPuCOvTByxg5CkZqXVx7bN8/vmRhx5ZeKI\nHYQktSryXucmjthBSFI7IlYfJo7YQUhSJyLM+6h84nDmuKTC69K8jxrOHLfikFQuXZr3UfmKo+jf\nQZLqdKH6MHHEDkKSQgg478PEETsISQol0Iq7Jo7YQUhSaBnP+zBxxA5Ckrohw+rDxBE7CEnqpgzm\nfZg4YgchSd3W4cgrE0fsICQpljbnfZg4YgchSTG1UX2YOGIHIUl50EL1UerEcRjJGlRXAndP8L6J\nQ5LGaVhxd7LdBkudOFYBI8B2TBwqnhouvqlIppl13tG1c2bn4U1rGHgO2Nrw+tkkCeFJYOkEn5sP\nbANi7ssrdaIWOwBV1723r7n1sKdGDuh9ZOu3tizsv2xgxYpd/ectPyp2XM06CziF/RPHFpLkMQvY\nARwDLAFWA68C/jz9/R7gi0ycHUfDhBxFrSRtdnrOdj/fyueaPXa646Z7/8om28m7Wona7fSc7Xy+\nlc80e2wzx710zMDCoSWnrb5q75zhtaODi4aupsNrZzcqjgeAxpmNPenPrwNPA18B3gLcAnwYeBb4\nk/T324AbKVeSmEitJG12es52P9/K55o9drrjWmmzyGolarfTc7bz+VY+0+yxzRz30jFfveP6Ww57\nauSAWY9s3bFlYf9lLcQzoW71D/QCdwEnp88HgYuBC9LnlwDHAVe0eN6dwIkZxCdJVfId4LXtfviA\nDAOJoe0vLklqTzduVU3kYWDuuOevBx6KFIskqQWxEsfu9OfZJLex5gPfiBSLJCln1pF0dr8APANc\nlL7eRzIcdyewLE5okiRJkiRJUigHA9cAnwbeHjkWqdEJwFrgjtiBSA3eSTJfbhh4c+RYuq6ffXND\nbowZiDQFE4fy6pUk//GeUqxRVa1oZa2rk0kmtgAc2pXoVHXtrsUmhdbO3+ZHgc+EDy28Vta6Ogd4\nT/p+Kb68cq+Vv88xVhzqhmb/No8mWUXkk8BAMycuQsXRylpXG4HTgDXA+m4FqEpr5e/zFcANwDyS\n/9lJITX7t3kGMERyq//dwAenO3FRlxw5nSRTjtlG8uXvBi6PEpG0z1R/n5dEiUhKTPa3eQXJf7ib\nUoSKQ5KUI0VNHK51pTzz71N5lcnfZlETh2tdKc/8+1ReVeZv07WulGf+fSqv/NuUJEmSJEmSJEmS\nJEmSJEmSJEkquBmxA5Byai/w+Ljn60iWnZYkaUIjAc5Z1NWopTpFXatKiuUp4GMk1ciXSfYRBzgE\nuAy4n2T59Fr6+oUkGzdtAO4BDgY+ATwBfA54EDiVZDmI1ePaeT9wTagvIUnK3v+R7JQ29liYvv49\nYFX6+8dJkgAkCWJs3Z9fZt/CcRcCu0gWlCM9zz8CB5HsWPki8CbgMJK1g16WHvcgycqlUu5YOksT\n+xnJtpsTuTn9eS/7Ese7SJLD2EJyRwGvGXfcU+nvC0gqjZ8DXyPZhQ3gp+lxv02y0c6BJFWJlDsm\nDql1Y9tx/oLkFhUkt33/iGRLzvHOAr7f5HnXklQx24HhDmOUgrGPQ8rGbSR7Nb88fT5WrTSOXLwH\nOJ/kVlUfMGvce98Efg1YRDKKS8olKw5pYoeS9G2M+WfgjxuOGU0fAJ8HjiFJDEcA3wXe0XAMJB3q\nbwAeI+lgf4J9t6sAbgd+nX0b7kiSKm4mycgqgNNJOsHH+xfgzK5GJEnKtcNJ9n3eCnwBOCN9/Uhg\nG/A3keKSJEmSJEmSJEmSJEmSJElt+X9xOLN5DU++cwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x4c58b10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "energies = np.logspace(0,2,100)\n",
    "\n",
    "for f in [lambda x:plain_numpy_powerlaw(x, 1, -2), sherpa_powerlaw, astromodels_powerlaw, astropy_powerlaw]:\n",
    "    \n",
    "    this_values = f(energies)\n",
    "    plt.plot(energies, this_values)\n",
    "\n",
    "plt.ylabel(\"Differential flux\")\n",
    "plt.xlabel(\"Energy\")\n",
    "plt.loglog()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Simple speed test</h1>\n",
    "\n",
    "Let's check how fast the 3 packages can evaluate the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 4.58 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "100000 loops, best of 10: 11.1 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r10 plain_numpy_powerlaw(energies, 1, -2.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 4.44 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "10000 loops, best of 10: 23.6 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r10 sherpa_powerlaw(energies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 10: 21 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r10 astromodels_powerlaw(energies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10000 loops, best of 10: 130 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r10 astropy_powerlaw(energies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sherpa and astromodels are very similar, and they are around 2x slower than the reference implementation. This is the price of all the nuts and bolts needed for a real model (access the current values of parameters, check for proper inputs and so on). Astropy however is 6x slower."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Simulations of a fit</h1>\n",
    "\n",
    "A better benchmark is to see how they behave in a more realistic situation. During a fit, the value of the parameters is changed a couple of times, and the function re-evaluated at each time. Let's build a function which \"simulates\" this situation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_iterations = 100\n",
    "\n",
    "def reference_fit():\n",
    "    \n",
    "    _ = plain_numpy_powerlaw(energies, 1.0, -2.0)\n",
    "\n",
    "def sherpa_fake_fit():\n",
    "    \n",
    "    # Iterations for this fake fit\n",
    "        \n",
    "    for i in range(n_iterations):\n",
    "        \n",
    "        # We use always the same values, to avoid the overhead of generating\n",
    "        # a random number. This will test anyway the speed of assigning to\n",
    "        # parameters\n",
    "        # Also, we assign directly to the .val attribute to avoid the overhead of the\n",
    "        # machinery which allows you to do sherpa_powerlaw.index = -2.0\n",
    "        \n",
    "        sherpa_powerlaw.index.val = -2.0\n",
    "        sherpa_powerlaw.ampl.val = 1.0\n",
    "        \n",
    "        _ = sherpa_powerlaw(energies)\n",
    "\n",
    "def astromodels_fake_fit():\n",
    "    \n",
    "    for i in range(n_iterations):\n",
    "        \n",
    "        # We use always the same values, to avoid the overhead of generating\n",
    "        # a random number. This will test anyway the speed of assigning to\n",
    "        # parameters\n",
    "        # Also, we assign directly to the .value attribute to avoid the overhead of the\n",
    "        # machinery which allows you to do astromodels_powerlaw.index = -2.0\n",
    "        \n",
    "        #astromodels_powerlaw.index.value = -2.0\n",
    "        #astromodels_powerlaw.logK.value = 0.0\n",
    "        \n",
    "        _ = astromodels_powerlaw(energies)\n",
    "\n",
    "def astropy_fake_fit():\n",
    "    \n",
    "    for i in range(n_iterations):\n",
    "        \n",
    "        # We use always the same values, to avoid the overhead of generating\n",
    "        # a random number. This will test anyway the speed of assigning to\n",
    "        # parameters\n",
    "        # Also, we assign directly to the .value attribute to avoid the overhead of the\n",
    "        # machinery which allows you to do astropy_powerlaw.alpha = 2.0\n",
    "        \n",
    "        astropy_powerlaw.alpha.value = 2.0\n",
    "        astropy_powerlaw.amplitude.value = 1.0\n",
    "        \n",
    "        _ = astropy_powerlaw(energies)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 8.26 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "100000 loops, best of 10: 10.4 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r10 reference_fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 10: 5.77 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r10 sherpa_fake_fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 10: 2.07 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r10 astromodels_fake_fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 loops, best of 10: 18.9 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit -r10 astropy_fake_fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here astromodels is 2.8x faster than sherpa, while astropy is 9x slower than astromodels, and 3.3x slower than sherpa. With respect to the reference implementation, astromodels is 19x slower, which is quite a lot. This is probably easy to reduce by using Cython, since almost all the overhead which generates that is given by the lookup of attributes and methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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": 0
}