2014-05-04 22:57:02 +00:00
{
"metadata": {
"name": "",
2014-05-05 19:34:22 +00:00
"signature": "sha256:3b47add80317bb8ad369eacd5528c941bb1a2cbd92ae5b33aae66d56ae35c741"
2014-05-04 22:57:02 +00:00
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
2014-05-05 00:16:02 +00:00
"[Sebastian Raschka](http://www.sebastianraschka.com) \n",
2014-05-04 22:57:02 +00:00
"last updated: 05/04/2014\n",
"\n",
2014-05-04 23:11:16 +00:00
"- [Link to this IPython Notebook on GitHub](https://github.com/rasbt/python_reference/blob/master/benchmarks/cython_least_squares.ipynb) \n",
2014-05-04 22:57:02 +00:00
"- [Link to the GitHub repository](https://github.com/rasbt/python_reference)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2014-05-05 19:34:22 +00:00
"#### The code in this notebook was executed in Python 3.4.0"
2014-05-04 22:57:02 +00:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<hr>\n",
"I am really looking forward to your comments and suggestions to improve and \n",
"extend this little collection! Just send me a quick note \n",
"via Twitter: [@rasbt](https://twitter.com/rasbt) \n",
"or Email: [bluewoodtree@gmail.com](mailto:bluewoodtree@gmail.com)\n",
"<hr>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Implementing the least squares fit method for linear regression and speeding it up via Cython"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name=\"sections\"></a>\n",
"<br>\n",
"<br>\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#Sections"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- [Introduction](#introduction)\n",
"- [Least squares fit implementations](#implementations)\n",
"- [Generating sample data and benchmarking](#sample_data)\n",
"- [Compiling the Python code via Cython in the IPython notebook](#cython_nb)\n",
2014-05-05 05:44:10 +00:00
"- [Performance growth rates for different sample sizes](#sample_sizes)\n",
2014-05-05 19:34:22 +00:00
"- [Bonus: How to use Cython without the IPython magic](#cython_bonus)\n",
"- [Appendix I: Cython vs. Numba](#numba)\n",
"- [Appendix II: Cython with and without type declarations](#type_declarations)"
2014-05-04 22:57:02 +00:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name=\"introduction\"></a>\n",
"<br>\n",
"<br>"
]
},
2014-05-05 19:34:22 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
" \n",
"(Note that this chart just reflects my rather objective thoughts after experimenting with Cython, and it is not based on real numbers or benchmarks.)\n",
"<br>\n",
"<br>\n",
"<br>\n",
"<br>"
]
},
2014-05-04 22:57:02 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[[back to top](#sections)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Linear regression via the least squares method is the simplest approach to performing a regression analysis of a dependent and a explanatory variable. The objective is to find the best-fitting straight line through a set of points that minimizes the sum of the squared offsets from the line. \n",
"The offsets come in 2 different flavors: perpendicular and vertical - with respect to the line. \n",
" \n",
" \n",
"\n",
2014-05-05 02:17:27 +00:00
"As Michael Burger summarizes it nicely in his article \"[Problems of Linear Least Square Regression - And Approaches to Handle Them](http://www.arsa-conf.com/archive/?vid=1&aid=2&kid=60101-220)\": \"the perpendicular offset method delivers a more precise result but is are more complicated to handle. Therefore normally the vertical offsets are used.\" \n",
"Here, we will also use the method of computing the vertical offsets.\n",
"\n"
2014-05-04 22:57:02 +00:00
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In more mathematical terms, our goal is to compute the best fit to *n* points $(x_i, y_i)$ with $i=1,2,...n,$ via linear equation of the form \n",
"$f(x) = a\\cdot x + b$. \n",
"Here, we assume that the y-component is functionally dependent on the x-component. \n",
"In a cartesian coordinate system, $b$ is the intercept of the straight line with the y-axis, and $a$ is the slope of this line."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to obtain the parameters for the linear regression line for a set of multiple points, we can re-write the problem as matrix equation \n",
"$\\pmb X \\; \\pmb a = \\pmb y$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\Rightarrow\\Bigg[ \\begin{array}{cc}\n",
"x_1 & 1 \\\\\n",
"... & 1 \\\\\n",
"x_n & 1 \\end{array} \\Bigg]$\n",
"$\\bigg[ \\begin{array}{c}\n",
"a \\\\\n",
"b \\end{array} \\bigg]$\n",
"$=\\Bigg[ \\begin{array}{c}\n",
"y_1 \\\\\n",
"... \\\\\n",
"y_n \\end{array} \\Bigg]$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With a little bit of calculus, we can rearrange the term in order to obtain the parameter vector $\\pmb a = [a\\;b]^T$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\Rightarrow \\pmb a = (\\pmb X^T \\; \\pmb X)^{-1} \\pmb X^T \\; \\pmb y$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"The more classic approach to obtain the slope parameter $a$ and y-axis intercept $b$ would be:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$a = \\frac{S_{x,y}}{\\sigma_{x}^{2}}\\quad$ (slope)\n",
"\n",
"\n",
"$b = \\bar{y} - a\\bar{x}\\quad$ (y-axis intercept)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where \n",
"\n",
"\n",
"$S_{xy} = \\sum_{i=1}^{n} (x_i - \\bar{x})(y_i - \\bar{y})\\quad$ (covariance)\n",
"\n",
"\n",
"$\\sigma{_x}^{2} = \\sum_{i=1}^{n} (x_i - \\bar{x})^2\\quad$ (variance)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name=\"implementations\"></a>\n",
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"## Least squares fit implementations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[[back to top](#sections)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1. The matrix approach"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, let us implement the equation:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\pmb a = (\\pmb X^T \\; \\pmb X)^{-1} \\pmb X^T \\; \\pmb y$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which I will refer to as the \"matrix approach\"."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"\n",
"def lin_lstsqr_mat(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" X = np.vstack([x, np.ones(len(x))]).T\n",
" return (np.linalg.inv(X.T.dot(X)).dot(X.T)).dot(y)"
],
"language": "python",
"metadata": {},
"outputs": [],
2014-05-05 19:34:22 +00:00
"prompt_number": 36
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2. The classic approach"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we will calculate the parameters separately, using only standard library functions in Python, which I will call the \"classic approach\"."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$a = \\frac{S_{x,y}}{\\sigma_{x}^{2}}\\quad$ (slope)\n",
"\n",
"\n",
"$b = \\bar{y} - a\\bar{x}\\quad$ (y-axis intercept)"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def classic_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" x_avg = sum(x)/len(x)\n",
" y_avg = sum(y)/len(y)\n",
" var_x = sum([(x_i - x_avg)**2 for x_i in x])\n",
" cov_xy = sum([(x_i - x_avg)*(y_i - y_avg) for x_i,y_i in zip(x,y)])\n",
" slope = cov_xy / var_x\n",
" y_interc = y_avg - slope*x_avg\n",
" return (slope, y_interc)"
],
"language": "python",
"metadata": {},
"outputs": [],
2014-05-05 19:34:22 +00:00
"prompt_number": 37
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3. Using the lstsq numpy function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For our convenience, `numpy` has a function that can also compute the leat squares solution of a linear matrix equation. For more information, please refer to the [documentation](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html)."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def numpy_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" X = np.vstack([x, np.ones(len(x))]).T\n",
" return np.linalg.lstsq(X,y)[0]"
],
"language": "python",
"metadata": {},
"outputs": [],
2014-05-05 19:34:22 +00:00
"prompt_number": 38
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4. Using the linregress scipy function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The last approach is using `scipy.stats.linregress()`, which returns a tuple of 5 different attributes, where the 1st value in the tuple is the slope, and the second value is the y-axis intercept, respectively. The documentation for this function can be found [here](http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.stats.linregress.html)."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import scipy.stats\n",
"\n",
"def scipy_lstsqr(x,y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" return scipy.stats.linregress(x, y)[0:2]"
],
"language": "python",
"metadata": {},
"outputs": [],
2014-05-05 19:34:22 +00:00
"prompt_number": 39
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name='sample_data'></a>\n",
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generating sample data and benchmarking"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[[back to top](#sections)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
2014-05-05 00:16:02 +00:00
"In order to test our different least squares fit implementation, we will generate some sample data: \n",
2014-05-04 22:57:02 +00:00
"- 500 sample points for the x-component within the range [0,500) \n",
"- 500 sample points for the y-component within the range [100,600) \n",
"\n",
"where each sample point is multiplied by a random value within\n",
"the range [0.8, 12)."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import random\n",
"random.seed(12345)\n",
"\n",
"x = [x_i*random.randrange(8,12)/10 for x_i in range(500)]\n",
"y = [y_i*random.randrange(8,12)/10 for y_i in range(100,600)]"
],
"language": "python",
"metadata": {},
"outputs": [],
2014-05-05 19:34:22 +00:00
"prompt_number": 40
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
2014-05-04 23:21:25 +00:00
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
2014-05-04 22:57:02 +00:00
"#### Visualization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To check how our dataset is distributed, and how the straight line, which we obtain via the least square fit method, we will plot it in a scatter plot. \n",
"Note that we are using our \"matrix approach\" here for simplicity, but after plotting the data, we will check whether all of the four different implementations yield the same parameters."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%pylab inline\n",
"from matplotlib import pyplot as plt\n",
"\n",
"slope, intercept = lin_lstsqr_mat(x, y)\n",
"\n",
"line_x = [round(min(x)) - 1, round(max(x)) + 1]\n",
"line_y = [slope*x_i + intercept for x_i in line_x]\n",
"\n",
"plt.figure(figsize=(8,8))\n",
"plt.scatter(x,y)\n",
"plt.plot(line_x, line_y, color='red', lw='2')\n",
"\n",
"plt.ylabel('y')\n",
"plt.xlabel('x')\n",
"plt.title('Linear regression via least squares fit')\n",
"\n",
"ftext = 'y = ax + b = {:.3f} + {:.3f}x'\\\n",
" .format(slope, intercept)\n",
"plt.figtext(.15,.8, ftext, fontsize=11, ha='left')\n",
"\n",
"plt.show()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAH4CAYAAACi3S9CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdUVMfbwPHv0lm6giAgooBRbNi7GCv2EnvXWOPPnhg1\nGksSxRhjoiYaE0SNsb4xEY0Fe+zGmijGLkZARRERFljYnfeP1Y0ELJRlAedzDkf2lpln7q48e++d\nO6MQQggkSZIkSSqSTIwdgCRJkiRJhiMTvSRJkiQVYTLRS5IkSVIRJhO9JEmSJBVhMtFLkiRJUhEm\nE70kSZIkFWEy0UtGd+jQIcqXL2/sMAqtSpUq8fvvv+drnSNHjuTTTz/N0b7e3t7s3bs3jyN6s/3y\nyy+UKlUKe3t7zp07Z5TPhFRwKeRz9FJ+8fb2JiQkhGbNmhk7FMmIypQpQ0hICE2bNjVI+QcOHKBf\nv378888/Bim/IPLx8eGrr76iffv2mdbNnDmT69ev8+OPPxohMqkgkGf0Ur5RKBQoFApjh6Gn0Wjy\nZJvXJYRAfq+Wnsmrz5YQgtu3b+Pv758n5UlFj0z0ktEdOHCAUqVK6V97e3uzYMECqlatiqOjIz17\n9iQ1NVW/ftu2bQQEBODk5ESDBg3466+/9OuCg4Px9fXF3t6eihUr8uuvv+rXrVy5kgYNGjBhwgSc\nnZ2ZNWtWplhmzpxJ165d6devHw4ODqxatYrHjx/z7rvv4u7ujqenJ9OnT0er1QKg1WqZOHEiLi4u\nlC1bliVLlmBiYqJf36RJE6ZNm0aDBg2wsbHh5s2b/P3337Ro0YLixYtTvnx5Nm3apK9/+/btVKxY\nEXt7ezw9PVmwYAEADx48oF27djg5OVG8eHEaN26c4Xg9uxSemprKuHHj8PDwwMPDg/Hjx6NWq/XH\n2dPTky+//BJXV1fc3d1ZuXJllu/Jhg0bqFWrVoZlCxcupGPHjgAMHDiQ6dOnA/Do0SPatWtHiRIl\nKFasGO3btycqKirLcv9LCKF/z5ydnenRowePHj3Sr+/WrRslS5bE0dGRwMBAIiIiXnisvvzyS1Qq\nFa1btyY6Oho7Ozvs7e25e/dupnpfdJwB5s+fr3+vV6xYgYmJCTdu3AB072dISIh+25UrV9KoUSP9\n67Fjx+Ll5YWDgwM1a9bk8OHD+nXZ/Wxdu3aNwMBAHB0dcXFxoWfPnpnakZqaip2dHRqNhqpVq+Ln\n5wf8+5nYuXMnc+fOZcOGDdjZ2VGtWrXXel+kIkZIUj7x9vYWe/fuzbR8//79wtPTM8N2derUETEx\nMSIuLk5UqFBBLFu2TAghxJkzZ0SJEiXEyZMnhVarFatWrRLe3t5CrVYLIYTYtGmTiImJEUIIsWHD\nBmFjYyPu3r0rhBAiNDRUmJmZiSVLlgiNRiOSk5MzxTJjxgxhbm4utmzZIoQQIjk5WXTq1EmMGDFC\nqFQqcf/+fVG7dm3x3XffCSGEWLp0qfD39xdRUVHi0aNHolmzZsLExERoNBohhBCBgYGidOnSIiIi\nQmg0GhEfHy88PT3FypUrhUajEWfPnhXOzs7i0qVLQggh3NzcxOHDh4UQQsTHx4szZ84IIYSYPHmy\nGDFihEhPTxfp6en6bf57XKdPny7q1asnYmNjRWxsrKhfv76YPn26/jibmZmJGTNmiPT0dLF9+3ah\nVCpFfHx8puOgUqmEnZ2duHr1qn5ZzZo1xYYNG4QQQgwcOFBf7sOHD8XmzZtFcnKyePLkiejWrZvo\n1KlT1h+C/8T71VdfiXr16omoqCihVqvF8OHDRa9evfTbhoaGisTERKFWq8W4ceNEQECAft2LjtWB\nAwcyfJ6y8qJ9d+zYIVxdXcXFixdFUlKS6NWrl1AoFOL69etCCCGaNGkiQkJCMsTXsGFD/es1a9aI\nuLg4odFoxIIFC4Sbm5tITU0VQmT/s9WzZ08xZ84cIYQQqamp4siRIy9sz/Mx/vcYz5w5U/Tr1++l\nx0Mq2uQZvVQgjRkzBjc3N5ycnGjfvj3nzp0DYPny5QwfPpxatWqhUCjo378/lpaWHDt2DICuXbvi\n5uYGQPfu3fHz8+PEiRP6ct3d3Rk1ahQmJiZYWVllWXf9+vXp0KEDAI8fP2bHjh0sXLgQa2trXFxc\nGDduHOvXrwdg48aNjBs3Dnd3dxwdHZkyZUqGy/MKhYKBAwdSoUIFTExM2LlzJ2XKlGHAgAGYmJgQ\nEBBAly5d2LhxIwAWFhZcvHiRhIQEHBwc9GdgFhYWxMTEcOvWLUxNTWnQoEGWsa9du5aPP/4YZ2dn\nnJ2dmTFjRoZ7s+bm5nz88ceYmprSunVrbG1tuXz5cqZyrK2t6dixI+vWrQPg6tWrXL58WX9cAH07\nixUrRufOnbGyssLW1papU6dy8ODBrN/Y//juu+/49NNPcXd3x9zcnBkzZvB///d/+rPagQMHYmNj\no193/vx5njx58tJjJV7j9siL9t24cSODBw/G398fpVKZ5VWfl+nTpw9OTk6YmJgwYcIEUlNTMxzf\n7Hy2LCwsuHXrFlFRUVhYWFC/fv1sxfKMkLeM3ngy0UsF0rNkDbqkk5iYCEBkZCQLFizAyclJ/3Pn\nzh1iYmIAWL16NdWqVdOvu3DhAg8fPtSX9fwtghfx9PTU/x4ZGUlaWholS5bUlzlixAhiY2MBiImJ\nyVDm8/tmVWdkZCQnTpzIEP/atWu5d+8eAD///DPbt2/H29ubJk2acPz4cQA++OADfH19admyJT4+\nPsybNy/L2KOjoyldurT+tZeXF9HR0frXxYsXx8Tk3//2SqVSf2z/q3fv3vpEv3btWn0y/y+VSsXw\n4cPx9vbGwcGBwMBAHj9+/FrJ5datW3Tu3Fl/LPz9/TEzM+PevXtoNBomT56Mr68vDg4OlClTBoVC\nwYMHD156rF7Hi/b97/vp5eX12mUCfPHFF/j7++Po6IiTkxOPHz/WxwvZ+2x9/vnnCCGoXbs2lSpV\nIjQ0NFuxSNIzZsYOQJJex7NOfF5eXnz00UdMnTo10zaRkZEMGzaMffv2Ua9ePRQKBdWqVct0hv2q\nep7fplSpUlhaWvLw4cMMCfKZkiVLZujdnVVP7+fL8/LyIjAwkPDw8Czrr1mzJr/++isajYbFixfT\nvXt3bt++ja2tLV988QVffPEFFy9epGnTptSuXZu33347w/7u7u7cunWLChUqAHD79m3c3d1f2uYX\nad68ObGxsZw/f57169fz1VdfZdmuBQsWcOXKFU6ePEmJEiU4d+4c1atXRwjxyuPt5eVFaGgo9erV\ny7Tuxx9/JCwsjL1791K6dGni4+MpVqyY/v180bF6nQ6fL9q3ZMmS3L59W7/d878D2NjYkJSUpH/9\n/P3/Q4cOMX/+fPbt20fFihUBMsT7/DGDV3+2XF1dWb58OQBHjhyhefPmBAYGUrZs2Ve273kFqQOs\nZBzyjF7KV2q1mpSUFP3P6/Y8fvbHcujQoSxbtoyTJ08ihCApKYnffvuNxMREkpKSUCgUODs7o9Vq\nCQ0N5cKFC9mK779noSVLlqRly5ZMmDCBJ0+eoNVquX79uv4Z5e7du/P1118THR1NfHw88+bNy/SH\n9fky27Vrx5UrV1izZg1paWmkpaXxxx9/8Pfff5OWlsZPP/3E48ePMTU1xc7ODlNTU0DXAfHatWsI\nIbC3t8fU1DTL5NCrVy8+/fRTHjx4wIMHD5g9ezb9+vXL1jF4xtzcnG7duvH+++/z6NEjWrRokaFN\nz9qVmJiItbU1Dg4OxMXFZety94gRI5g6dao+ocbGxhIWFqYv19LSkmLFipGUlJThy93LjpWrqysP\nHz4kISEhyzpftm/37t1ZuXIlly5dQqVSZWpLQEAAmzdvJjk5mWvXrhESEqJ/v588eYKZmRnOzs6o\n1Wpmz579whjg1Z+tTZs2cefOHQAcHR1RKBRZvuev4ubmxq1bt+Tl+zeYTPRSvmrTpg1KpVL/M2vW\nrFc+dvf8+ho1avD999/zv//9
"text": [
2014-05-05 19:34:22 +00:00
"<matplotlib.figure.Figure at 0x10f9e8a50>"
2014-05-04 22:57:02 +00:00
]
}
],
2014-05-05 19:34:22 +00:00
"prompt_number": 42
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
2014-05-04 23:21:25 +00:00
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
2014-05-04 22:57:02 +00:00
"#### Comparing the results from the different implementations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As mentioned above, let us confirm that the different implementation computed the same parameters (i.e., slope and y-axis intercept) as solution for the linear equation."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import prettytable\n",
"\n",
"params = [appr(x,y) for appr in [lin_lstsqr_mat, classic_lstsqr, numpy_lstsqr, scipy_lstsqr]]\n",
"\n",
"print(params)\n",
"\n",
"fit_table = prettytable.PrettyTable([\"\", \"slope\", \"y-intercept\"])\n",
"fit_table.add_row([\"matrix approach\", params[0][0], params[0][1]])\n",
"fit_table.add_row([\"classic approach\", params[1][0], params[1][1]])\n",
"fit_table.add_row([\"numpy function\", params[2][0], params[2][1]])\n",
"fit_table.add_row([\"scipy function\", params[3][0], params[3][1]])\n",
"\n",
"print(fit_table)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[array([ 0.95181895, 107.01399744]), (0.95181895319126741, 107.01399744459181), array([ 0.95181895, 107.01399744]), (0.95181895319126764, 107.01399744459175)]\n",
"+------------------+----------------+---------------+\n",
"| | slope | y-intercept |\n",
"+------------------+----------------+---------------+\n",
"| matrix approach | 0.951818953191 | 107.013997445 |\n",
"| classic approach | 0.951818953191 | 107.013997445 |\n",
"| numpy function | 0.951818953191 | 107.013997445 |\n",
"| scipy function | 0.951818953191 | 107.013997445 |\n",
"+------------------+----------------+---------------+\n"
]
}
],
2014-05-05 19:34:22 +00:00
"prompt_number": 43
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
2014-05-04 23:21:25 +00:00
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
2014-05-04 22:57:02 +00:00
"#### Initial performance comparison"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For a first impression how the performances of the different least squares implementations compare against each other, let us do a quick benchmark using the `timeit` module via IPython's `%timeit` magic."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import timeit\n",
"\n",
"for lab,appr in zip([\"matrix approach\",\"classic approach\",\n",
" \"numpy function\",\"scipy function\"],\n",
" [lin_lstsqr_mat, classic_lstsqr, \n",
" numpy_lstsqr, scipy_lstsqr]):\n",
" print(\"\\n{}: \".format(lab), end=\"\")\n",
" %timeit appr(x, y)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
2014-05-05 19:34:22 +00:00
"matrix approach: 10000 loops, best of 3: 163 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"\n",
2014-05-05 19:34:22 +00:00
"classic approach: 1000 loops, best of 3: 1.55 ms per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"\n",
2014-05-05 19:34:22 +00:00
"numpy function: 1000 loops, best of 3: 221 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"\n",
2014-05-05 19:34:22 +00:00
"scipy function: 1000 loops, best of 3: 362 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
}
],
2014-05-05 19:34:22 +00:00
"prompt_number": 44
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"The timing above indicates, that the \"classic\" approach (Python's standard library functions only) is significantly worse in performance than the other implemenations - roughly by a magnitude of 10."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>\n",
"<a name=\"cython_nb\"></a>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compiling the Python code via Cython in the IPython notebook"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[[back to top](#sections)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Maybe we can speed things up a little bit via [Cython's C-extensions for Python](http://cython.org). Cython is basically a hybrid between C and Python and can be pictured as compiled Python code with type declarations. \n",
"Since we are working in an IPython notebook here, we can make use of the IPython magic: It will automatically convert it to C code, compile it, and load the function. \n",
2014-05-05 19:34:22 +00:00
"Just to be thorough, let us also compile the other functions, which already use numpy objects.\n",
"\n",
"**Note** \n",
"Of course Cython has much more horsepower under its hood - more than I am showing in this article (for example, I am not using Cython's type definitions via `cdef` here). Here, I want to focus on how to speed up existing Python code by making only minimal changes to it. \n",
"[In a later section - Appendix II](#type_declarations) We will see how static type declarations can further improve the performance via Cython."
2014-05-04 22:57:02 +00:00
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%load_ext cythonmagic"
],
"language": "python",
"metadata": {},
2014-05-05 19:34:22 +00:00
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"The cythonmagic extension is already loaded. To reload it, use:\n",
" %reload_ext cythonmagic\n"
]
}
2014-05-05 15:40:43 +00:00
],
2014-05-05 19:34:22 +00:00
"prompt_number": 45
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%cython\n",
"import numpy as np\n",
"import scipy.stats\n",
"cimport numpy as np\n",
"\n",
"def cy_lin_lstsqr_mat(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" X = np.vstack([x, np.ones(len(x))]).T\n",
" return (np.linalg.inv(X.T.dot(X)).dot(X.T)).dot(y)\n",
"\n",
"def cy_classic_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" x_avg = sum(x)/len(x)\n",
" y_avg = sum(y)/len(y)\n",
" var_x = sum([(x_i - x_avg)**2 for x_i in x])\n",
" cov_xy = sum([(x_i - x_avg)*(y_i - y_avg) for x_i,y_i in zip(x,y)])\n",
" slope = cov_xy / var_x\n",
" y_interc = y_avg - slope*x_avg\n",
" return (slope, y_interc)\n",
"\n",
"def cy_numpy_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" X = np.vstack([x, np.ones(len(x))]).T\n",
" return np.linalg.lstsq(X,y)[0]\n",
"\n",
"def cy_scipy_lstsqr(x,y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" return scipy.stats.linregress(x, y)[0:2]"
],
"language": "python",
"metadata": {},
"outputs": [],
2014-05-05 19:34:22 +00:00
"prompt_number": 46
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Comparing the compiled Cython code to the original Python code"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import timeit\n",
"\n",
"for lab,appr in zip([\"matrix approach\",\"classic approach\",\n",
" \"numpy function\",\"scipy function\"],\n",
" [(lin_lstsqr_mat, cy_lin_lstsqr_mat), \n",
" (classic_lstsqr, cy_classic_lstsqr),\n",
" (numpy_lstsqr, cy_numpy_lstsqr),\n",
" (scipy_lstsqr, cy_scipy_lstsqr)]):\n",
" print(\"\\n\\n{}: \".format(lab))\n",
" %timeit appr[0](x, y)\n",
" %timeit appr[1](x, y)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"\n",
"matrix approach: \n",
2014-05-05 19:34:22 +00:00
"10000 loops, best of 3: 165 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
2014-05-05 19:34:22 +00:00
"10000 loops, best of 3: 166 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"\n",
"\n",
"classic approach: \n",
2014-05-05 19:34:22 +00:00
"1000 loops, best of 3: 1.59 ms per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
2014-05-05 19:34:22 +00:00
"10000 loops, best of 3: 127 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"\n",
"\n",
"numpy function: \n",
2014-05-05 19:34:22 +00:00
"1000 loops, best of 3: 220 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
2014-05-05 19:34:22 +00:00
"1000 loops, best of 3: 221 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"\n",
"\n",
"scipy function: \n",
2014-05-05 19:34:22 +00:00
"1000 loops, best of 3: 361 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
2014-05-05 19:34:22 +00:00
"1000 loops, best of 3: 367 \u00b5s per loop"
2014-05-04 22:57:02 +00:00
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
}
],
2014-05-05 19:34:22 +00:00
"prompt_number": 47
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>\n",
"As we've seen before, our \"classic\" implementation of the least square method is pretty slow compared to using numpy's functions. This is not surprising, since numpy is highly optmized and uses compiled C/C++ and Fortran code already. This explains why there is no significant difference if we used Cython to compile the numpy-objects-containing functions. \n",
"However, we were able to speed up the \"classic approach\" quite significantly, roughly by 1500%.\n",
"\n",
"The following 2 code blocks are just to visualize our results in a bar plot."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import timeit\n",
"\n",
"funcs = ['classic_lstsqr', 'cy_classic_lstsqr', \n",
" 'lin_lstsqr_mat', 'numpy_lstsqr', 'scipy_lstsqr']\n",
"labels = ['classic approach','classic approach (cython)', \n",
" 'matrix approach', 'numpy function', 'scipy function']\n",
"\n",
"times = [timeit.Timer('%s(x,y)' %f, \n",
" 'from __main__ import %s, x, y' %f).timeit(1000)\n",
" for f in funcs]\n",
"\n",
"times_rel = [times[0]/times[i+1] for i in range(len(times[1:]))]"
],
"language": "python",
"metadata": {},
"outputs": [],
2014-05-05 19:34:22 +00:00
"prompt_number": 50
2014-05-04 22:57:02 +00:00
},
{
"cell_type": "code",
"collapsed": false,
"input": [
2014-05-05 19:34:22 +00:00
"#%pylab inline\n",
"#import matplotlib.pyplot as plt\n",
2014-05-04 22:57:02 +00:00
"\n",
"x_pos = np.arange(len(funcs))\n",
"plt.bar(x_pos, times, align='center', alpha=0.5)\n",
"plt.xticks(x_pos, labels, rotation=45)\n",
"plt.ylabel('time in ms')\n",
"plt.title('Performance of different least square fit implementations')\n",
2014-05-05 19:34:22 +00:00
"plt.grid()\n",
2014-05-04 22:57:02 +00:00
"plt.show()\n",
"\n",
"x_pos = np.arange(len(funcs[1:]))\n",
"plt.bar(x_pos, times_rel, align='center', alpha=0.5, color=\"green\")\n",
"plt.xticks(x_pos, labels[1:], rotation=45)\n",
"plt.ylabel('relative performance gain')\n",
"plt.title('Performance gain compared to the classic least square implementation')\n",
2014-05-05 19:34:22 +00:00
"plt.grid()\n",
2014-05-04 22:57:02 +00:00
"plt.show()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
2014-05-05 19:34:22 +00:00
"png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAFhCAYAAABwNN3iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlcFPX/B/DXIiCgHAJeHIIIJiKCigeagmkqKGhqKnnh\nkWhqmh2alqJ5YZ5FmWnKN1M0yxMBTW0xD8Ajj0TwalFRSLwQQWGX9+8P2vmxHALmsruz7+fj4UNm\nZ3bm/Z6ZnffO5zMzKyEiAmOMMb1noOkAGGOMaQcuCIwxxgBwQWCMMfYvLgiMMcYAcEFgjDH2Ly4I\njDHGAIi8IGRlZaFbt26wsLDAxx9/rOlwNC4/Px9BQUGwsrLC0KFDK51eKpXC0dFRGG7VqhWOHj0K\nACAijBkzBtbW1ujUqRMAYO3atWjYsCEsLCzw8OFD9STxkkrnwqouLS0N3t7esLCwwNdff41JkyZh\n4cKFVX5/daevDmdnZxw+fFgt89akLVu2oHfv3jW/YNIyTk5OZGpqSnXr1qWGDRtSaGgo5ebmvtS8\nFixYQIMGDXrFEequH3/8kTp06EAKhaJK0//+++/k4OBQ7rijR4+Sg4MD5eXlERFRQUEBmZqa0sWL\nF19ZvNXh5OREhw8frnD8i3J5VUaPHk2fffaZWpehCWPHjqUZM2aUO64m1uuLODs7v3C7a8KmTZvo\n9ddfr/L0f//9N0kkkip/LtVJ684QJBIJYmJi8OTJE5w9exanT5+u9rcLIkJRURHS09Ph7u7+UnHI\n5fKXep82S09PR/PmzWFg8N83e3p6OpydnWFqagoAyMzMxLNnz156fRcVFf2neCQSCYjvsazUy+zX\n6enpaNmypRqiYSVpxf6r4YJURumK/9FHH1G/fv2IiOjkyZPk6+tLVlZW5OXlRVKpVJjOz8+P5syZ\nQ126dCFTU1MaMWIEGRkZkbGxMdWtW5cOHz5Mz58/p2nTppGdnR3Z2dnR9OnT6fnz50RU/E3H3t6e\nIiIiqFGjRjRy5EgKDw+nwYMH04gRI8jc3Jw8PT3pypUrtHjxYmrQoAE1adKEDh48KMSwceNGcnd3\nJ3Nzc3JxcaF169YJ45TzX7FiBTVo0IAaN25MmzZtEsbn5eXRjBkzyMnJiSwtLen111+n/Pz8SvMu\nLSUlhfz8/MjKyoo8PDxo7969REQ0d+5cMjY2JiMjI6pbty5t3LixzHvz8vJo9OjRVK9ePWrZsiUt\nW7ZM5dufk5MTHTp0iDZs2EAmJiZUq1Ytqlu3LoWEhFCdOnVIIpFQ3bp1qUePHkREdPnyZerZsydZ\nW1vTa6+9Rj///LMwr9GjR9PEiRMpICCA6tSpQ4cPH6aMjAwaOHAg1a9fn5o2bUpfffWVMP28efPo\n7bffplGjRpG5uTl5eHjQ6dOniYhoxIgRZGBgIJxZfvnll2VyK/1N9kXLSkpKok6dOpGVlRU1btyY\npkyZQgUFBcL46dOnU4MGDcjCwoI8PT3pr7/+onXr1qnsb8HBweVun/LeS0SUnZ1NQUFBZGFhQR06\ndKDPPvtM+JZZ3jdIPz8/2rBhAxERXbt2jbp37042NjZka2tLw4cPp0ePHqlst4iICPL09CQTExNS\nKBRV3qe6d+9OtWrVIhMTEzI3N6crV64IZ0JPnz4lExMTMjAwoLp165K5uTndvXu3zDxKnjkpPwfL\nli2j+vXrU+PGjWnXrl20f/9+cnNzI2tra1qyZInKdh80aBANHTqUzM3NqW3btnT+/HlhfMnjRVFR\nES1ZsoSaNWtGNjY2NGTIEHrw4IHKOty0aRM5OjqStbU1rV27lpKTk8nT05OsrKxoypQpKnH/8MMP\n5O7uTvXq1aPevXtTenq6ME4ikdB3331Hbm5uZGVlRZMnTyai4s9fyc9GvXr1iIgoJiaGvL29ycLC\nghwdHSk8PFyYl6Ojo/DZMTc3p5MnT5Y5yzh+/Dj5+PiQpaUltW/fnk6cOKGyL3z++efUpUsXMjc3\np169elF2djYREeXn59Pw4cPJxsaGrKysqH379pSVlVXutiYi0sqCcOjQISIiunnzJnl4eNDcuXPp\n9u3bZGNjQ3FxcURE9Ntvv5GNjY2QuJ+fHzk5OVFKSgopFAoqLCyk0NBQ+vzzz4V5f/755+Tr60v3\n7t2je/fuUefOnYXxv//+OxkaGtKsWbOooKCA8vPzad68eWRiYkIHDx4kuVxOo0aNIicnJ1q8eDHJ\n5XJav349NW3aVJj//v376caNG0RElJCQQGZmZnT27FmV+c+bN4/kcjnFxsaSmZmZ8MF97733qHv3\n7nTnzh3hA/v8+fMK8753716ZdVdQUEDNmjWjJUuWUGFhIR05coTMzc0pLS2NiIjCw8Np5MiRFa77\nmTNnUrdu3ejhw4d069Yt8vDwIEdHR5Vto/zwRUVFqeywMplM5aCVm5tLDg4OFBUVRQqFgv7880+y\ntbWllJQUIio+SFhaWgo7dl5eHrVt25a++OILKiwspBs3bpCLiwsdOHCAiEjYFnFxcVRUVESffvop\nderUqdzYylOyICgUihcu68yZM5SUlEQKhYJkMhm5u7vT6tWriYgoPj6e2rVrR48fPyYiotTUVOEg\nWHp/K+1F7x06dCgNHTqU8vLy6K+//iJ7e3vq2rUrEZVfEPz9/emHH34gouKCcOjQISooKKB79+5R\nt27daPr06cK0Tk5O1KZNG7p9+zY9e/asWvtU6WWVzlMqlVbaZFRyeuXn4IsvvhA+QzY2NvTOO+9Q\nbm4uXbp0iUxNTUkmkxFR8XY3MjKiX3/9leRyOS1fvpyaNm1KcrmciFS3++rVq8nX15cyMjKooKCA\nwsLCKCQkRGUdTpo0iZ4/f04HDx4kY2NjGjBgAN27d48yMjKoQYMGlJCQQEREu3fvJldXV0pNTSWF\nQkELFy6kzp07CzlJJBIKCgqix48f082bN6l+/foUHx9PRGU/G8r1pCz+Fy5coIYNG9Lu3buJqOxn\nh0i12en+/ftkZWVFP/30EykUCoqOjqZ69eoJxc7Pz49cXV3p6tWrlJ+fT/7+/jRr1iwiIvruu+8o\nKCiI8vPzqaioiM6ePUs5OTkVbiutKwhOTk5Ut25dsrKyIicnJ5o8eTLl5+fT0qVLyxzMevfuTf/7\n3/+IqHinnTdvnsr40NBQlTbdZs2aCR8CIqIDBw6Qs7MzERXvqMbGxsIZA1HxztirVy9heO/evVS3\nbl0qKioiIqKcnBySSCTCB7y0AQMG0Jo1a4T5m5qaqmz0Bg0aCAceU1NTunDhQpl5VJZ3SUePHqVG\njRqpvBYSEiJ8G5k3bx6NGDGi3FiJSOWgSET0/fffq3zYS374Sn+DKX3Q2rZtm3BAU5owYQLNnz+f\niIoLwujRo4VxiYmJ1KRJE5XpFy9eTGPGjBFif/PNN4VxygNHebGVp2RBqGxZpa1atYreeustIiI6\nfPgwNW/enBITE8u0+Zbe30o7cuRIue+Vy+VkZGQkFG4iotmzZ7/wDKH0QbqkXbt2UZs2bYRhZ2dn\nlbPR6uxTymUpz0ZK51mVPoTS05uampb5DCUnJwvTt2vXjvbs2UNExdvd19dXGFdUVESNGzemY8eO\nCbkpt7u7u7vKPnDnzh0yMjIihUIhrMM7d+4I421sbFTOWgcNGiR8Xvv06aOyfhUKBZmZmdHNmzeJ\nqLggHD9+XBg/ZMgQWrp0KRFVrQ9h2rRp9MEHHxBR+du35Dx+/PFH6tixo8r7fX19KSoqioiKt8+i\nRYuEcd9++y316dOHiIpbLTp37lzusaU8hppusipNIpFgz549eOONN1ReT09Px44dO7Bv3z7hNblc\nrjJdZVeR3LlzB05OTsJwkyZNcOfOHWG4fv36MDY2VnlPgwYNhL9NTU1ha2sLiUQiDANAbm4uLCws\nEBcXh/nz5+Pq1asoKipCXl4eWrduLbzfxsZGpf3ezMwMubm5yM7OxrNnz9CsWbMyMVcl75L5lV4H\nTk5OyMjIeOF6qej9TZo0qdL7
2014-05-04 22:57:02 +00:00
"text": [
2014-05-05 19:34:22 +00:00
"<matplotlib.figure.Figure at 0x10eeafd50>"
2014-05-04 22:57:02 +00:00
]
},
{
"metadata": {},
"output_type": "display_data",
2014-05-05 19:34:22 +00:00
"png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAFhCAYAAAALEB8uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYFGfXBvB7qNJ7EaSIxgiKghWxYVBj16BGiSWoiSXR\naEyR1xjljbFrLNGYmESJPTHGGuyK+tl7I4JtsYIFlSqwy/n+8N2JC4iKw8yynN91eSWz7Mw8c+/O\nnJ3nmdkViIjAGGOMGTgjpRvAGGOMyYELHmOMsQqBCx5jjLEKgQseY4yxCoELHmOMsQqBCx5jjLEK\nQdaCl5qaihYtWsDW1hZffPGFnKsulzp06IBly5Yp3Qy95evri127dsm2PiMjI1y9erVM1xEVFYWv\nv/66zJZvY2MDlUr1SvPEx8fDy8urbBpk4GrXro19+/ZJvlyVSgUjIyMUFBRIvmylDRs2DN9++22Z\nLNvkRU/w9fXF3bt3YWxsDCsrK7Rv3x7z58+HlZXVK69s0aJFcHV1RXp6eqkaW9HExcUp3QS9JggC\nBEEo9m9RUVHw8vLCxIkTS7XssLAw9OvXD4MGDXqdJr6ykrZJChkZGWW27Nfxuq+Xvjp//rzSTZDV\nq76OsbGx+PXXX7F//37xsYULF5ZV8158hicIAjZv3oyMjAycPHkSx48ff+XqS0QoKChAcnIy/P39\nS9VQtVpdqvlY+aBvr29ZFp0X4e+C0E/69h5lpUAv4OvrS7t27RKnP//8c+rUqRMRER06dIiaNGlC\n9vb2VLduXYqPjxef17JlS/rqq6+oadOmZGFhQX379iVTU1MyMzMja2tr2rVrF+Xm5tLIkSPJw8OD\nPDw8aNSoUZSbm0tERHv27CFPT0+aNm0aubu7U79+/SgmJoZ69OhBffv2JRsbGwoMDKSkpCSaPHky\nubq6kre3N23fvl1sw+LFi8nf359sbGzIz8+PfvrpJ/Fv2uXPmjWLXF1dqXLlyrRkyRLx79nZ2TR6\n9Gjy8fEhOzs7atasGeXk5Lxwuws7ceIEBQUFkY2NDfXs2ZPeffddGjduHBERpaWlUceOHcnFxYUc\nHByoU6dOdPPmTZ0Mf/nlFyIiWrJkCTVt2pQ+//xzcnBwoKpVq9KWLVueu97r16/TO++8Qy4uLuTk\n5ETDhw8nIiKNRkMTJ04kHx8fcnV1pf79+9Pjx4+JiOjatWskCAItWbKEvLy8yNHRkRYuXEhHjx6l\nwMBAsre3F5ejbVNoaCgNHz6c7OzsqGbNmjrvlZfJX/v69u/fnwoKCmjKlClUrVo1cnJyonfffZfS\n0tLEeZYuXUre3t7k5OREkyZNKvLe1Prpp5903mtdunQhIqKEhARq2bIl2dvbU61atWjjxo3FZjd2\n7FgyNjamSpUqkbW1NY0YMYKIiARBoB9//JHeeOMNsre3p48//lhnvl9//ZX8/f3JwcGB3n77bUpO\nTn7u67N//37xPeTl5UW//fYbERFFRUW99PtjyZIl5OfnRzY2NlS1alVasWIFERFdunSJWrRoQXZ2\nduTs7Ey9evUS5xEEga5cuUJEJb/Hn7Vnzx6qUqWKOH3r1i2KiIggFxcXqlq1Ks2bN0/825EjRygk\nJITs7e2pcuXKNHz4cMrLyxP/PmrUKHJ1dSVbW1sKDAyk8+fPP/f1Kqy4eYmI7t+/T507dyZbW1tq\n1KgRjRs3jpo1a0ZE/76nNRqNuJxn96vLly9Tq1atyMnJiZydnalPnz706NEj8bk+Pj40bdo0CgwM\npEqVKpFGo3ml/d/Hx0d8j06YMOGVjl8tW7ak6OhoatSoEdna2lLXrl3F/aHwdj169IgGDhxIlStX\nJk9PTxo3bpz4N+1++umnn5K9vT1Vq1aNDhw4QIsXLyYvLy9ydXUV339ERE+ePKHPPvuMvL29yc3N\njYYOHSq+L0o6bj7vddTu0zY2NhQQEEDr1q0joqf7Y6VKlcjY2Jisra3JwcGBiIjef/99cR8gIlq0\naBFVr16dHB0dqUuXLnT79m3xby/aJwt7qYK3c+dOInp6EK1VqxaNHz+ebt68SU5OTuJBd8eOHeTk\n5ET3798XXywfHx9KSEggjUZD+fn5FBUVRV9//bW47K+//pqaNGlC9+7do3v37lFoaKj49z179pCJ\niQlFR0dTXl4e5eTk0IQJE6hSpUq0fft2UqvV1L9/f/Lx8aHJkyeTWq2mn3/+mapWrSou/++//6ar\nV68SEdHevXvJ0tKSTp48qbP8CRMmkFqtpri4OLK0tBTf7B999BG1atWKbt++Lb7Jc3Nzn7vd9+7d\nK5Jdbm4ueXt707x580itVtNff/1FZmZm4jY+ePCA/vrrL8rJyaGMjAzq2bMndevWTZw/LCyMfv31\nVyJ6+qY1NTWlX375hQoKCmjhwoXk4eFR7GumVqupTp06NHr0aMrOzqYnT57QgQMHiOjpQbl69ep0\n7do1yszMpIiICOrXrx8R/bsTDRs2jHJzc2n79u1kZmZG3bp1o3v37tGtW7fI1dWV9u7dK7bJxMSE\n5syZQ2q1mn7//Xeys7MTd8qXyf/Z13fOnDnUpEkTunXrFuXl5dGQIUMoMjKSiIguXLhA1tbWtH//\nfsrNzaXRo0eTiYlJsQWPiIq81/Ly8qhatWo0ZcoUys/Pp927d5ONjQ0lJiYWO/+z2WsJgkCdO3em\nx48f0/Xr18nFxYW2bt1KRETr16+n6tWr08WLF0mj0dC3335LoaGhxS5bpVKRjY0NrV69mtRqNT14\n8IBOnz4ttlu7s5f0/sjMzCRbW1tKSkoiIqKUlBS6cOECERH17t2bJk+eTERP34Pa1167DdqC97z3\neGHPFjyNRkP16tWjiRMnUn5+Pl29epX8/Pxo27ZtRPT0A96RI0dIo9GQSqUif39/mjNnDhERbd26\nlerXry9+wLp48SLduXOn2NersJLm7dWrF/Xq1Yuys7Pp/Pnz5OnpSc2bNyei4gves6/t5cuXaefO\nnZSXl0f37t2jFi1a0KhRo8Tn+vj4UHBwMN28eZOePHnySvs/ke4Jw6sev1q2bEmenp504cIFysrK\nou7du1Pfvn2L3a5u3brR0KFDKTs7m+7evUuNGjUSP2Bq99PY2FgqKCigcePGkaenp/hhZPv27WRj\nY0NZWVlE9PSDRdeuXenhw4eUkZFBnTt3pv/85z/ie6Gk42Zxr+OaNWvE1+r3338nKysrSklJISKi\n2NhY8cOJ1rPL2LVrFzk7O9OpU6coNzeXRowYQS1atBCfW9I+WZwXFjwfHx+ytrYme3t78vHxoY8/\n/phycnJo6tSp4oFS6+233xY/KYSFhdGECROKbMizlbtatWo6Zynbtm0jX19fInoarJmZmc4OOGHC\nBGrbtq04vXHjRrK2tqaCggIiIkpPTydBEMSdorBu3brR3LlzxeVbWFjo7Aiurq7izmphYUFnz54t\nsowXbfez9u7dS56enjqPNWvW7Lk79qlTp8RPOURFC1716tXFv2VlZZEgCJSamlpkOQcPHiQXFxed\nbdN66623aOHCheJ0YmIimZqakkajEXeiZz9BOTk50R9//CFOd+/eXTyALVmypEjRbdSoES1btqzY\n7Sucf+HX19/fX6eA3b59m0xNTUmtVtN///tfsfhpt9/MzKzEgvfse23fvn3k7u6u85zIyEiKiYkp\ndv6wsDDxLEBLEASd4vHuu+/StGnTiIioXbt2OgVSo9GQpaUlXb9+vciyJ0+eTBERES/V7mc9+/7I\nzMwke3t7Wrt2LWVnZ+s8r3///jR48GCds8Fnt+HKlSslvscLe7bgHT58mLy9vYtsz4ABA4qdd/bs\n2fTOO+8Q0dODV40aNejw4cNF3pslbTcR0e7du4udV61Wk6mpqc4Hl7Fjx5Z4hlfchxmtdevWUXBw\nsDjt6+ur0/PzKvu/dv5nC96rHL/CwsLEQkP09IzIzMyMCgoKdLYrJSWFzM3Ndc7OV65cSa1atSKi\np/vpG2+8If7t7NmzJAgC3b17
2014-05-04 22:57:02 +00:00
"text": [
2014-05-05 19:34:22 +00:00
"<matplotlib.figure.Figure at 0x10eeac310>"
2014-05-04 22:57:02 +00:00
]
}
],
2014-05-05 19:34:22 +00:00
"prompt_number": 51
2014-05-04 22:57:02 +00:00
},
2014-05-05 05:44:10 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name=\"sample_sizes\"></a>\n",
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Performance growth rates for different sample sizes"
]
},
2014-05-05 19:34:22 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[[back to top](#sections)]"
]
},
2014-05-05 05:44:10 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the plot above, we've seen how the different implemantations perform for a fixed sample size n=500. Now, let us take a look at the effect of the sample size on the relative performances for each approach. We will consider the sample sizes 10, 100, 1000, 10000, 100000, and 1000000."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import timeit\n",
2014-05-05 19:34:22 +00:00
"import random\n",
"random.seed(12345)\n",
2014-05-05 05:44:10 +00:00
"\n",
"funcs = ['cy_classic_lstsqr', \n",
" 'lin_lstsqr_mat', 'numpy_lstsqr', 'scipy_lstsqr']\n",
"labels = ['classic approach (cython)', 'matrix approach', \n",
" 'numpy function', 'scipy function']\n",
"orders_n = [10**n for n in range(1, 7)]\n",
"times_n = {f:[] for f in funcs}\n",
"\n",
"for n in orders_n:\n",
" x = [x_i*random.randrange(8,12)/10 for x_i in range(n)]\n",
" y = [y_i*random.randrange(10,14)/10 for y_i in range(n)]\n",
" for f in funcs:\n",
" times_n[f].append(timeit.Timer('%s(x,y)' %f, \n",
" 'from __main__ import %s, x, y' %f).timeit(1000))"
],
"language": "python",
"metadata": {},
"outputs": [],
2014-05-05 19:34:22 +00:00
"prompt_number": 22
2014-05-05 05:44:10 +00:00
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#%pylab inline\n",
"#import matplotlib.pyplot as plt\n",
"\n",
"plt.figure(figsize=(10,8))\n",
"\n",
"for f in times_n.keys():\n",
" plt.plot(orders_n, times_n[f], alpha=0.5, label=f, marker='o', lw=2)\n",
"\n",
"plt.xlabel('sample size n')\n",
"plt.ylabel('time in ms')\n",
"plt.xlim([0,max(orders_n) + max(orders_n) * 0.1])\n",
"plt.legend(loc=2)\n",
2014-05-05 19:34:22 +00:00
"plt.grid()\n",
2014-05-05 05:44:10 +00:00
"\n",
"plt.title('Performance of least square fit implementations for different sample sizes')\n",
"plt.show()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
2014-05-05 19:34:22 +00:00
"png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAH4CAYAAAACdDpdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlcVOX+B/DPDLIz7JvsCl4VV9IQxQWFXCqXXFBLBDW9\nppZbtyxTKXOrvLef3bxXCyXF3VLLJU2QUivJzA0MFQUUUFAUZZEZZr6/P+YyMayDMMyZme/79fIl\nZ+ac5zznfM+c+c5znvMcERERGGOMMcZYixLrugKMMcYYY8aIkzDGGGOMMR3gJIwxxhhjTAc4CWOM\nMcYY0wFOwhhjjDHGdICTMMYYY4wxHeAkzIDdvXsX/fv3h62tLf7xj3/oujo6V1ZWhuHDh8Pe3h7j\nx4+v8X5sbCyioqJ0UDP9t2/fPnh7e8PW1hbnz59H586d8dNPP2m8fGPn11RmZibEYjEUCkWzl61r\nr732Gj788MMWXWdlnCUSCS5cuNDs5Vf9DGZnZ0MikaByFKXazmdTpkyBo6MjQkJCmr0uhiA5ORne\n3t7NWubJkyfRoUOHZi3TmLXSdQWYOj8/P+Tn58PExATW1tYYNmwY/v3vf8Pa2rrRZW3cuBGurq54\n9OiRFmqqf/bu3Yv8/HwUFhZCLK75+0MkEmm9Dn5+fti0aRMGDRqk9XW1pDfffBPr16/H8OHDAQCX\nL19WvRcbG4uMjAxs3bq1zuWrzm8MYmJi4O3tjeXLl2s0f3x8POLi4nDy5EnVa//5z3+0Vb06VY9z\nc6v6GfTx8cHjx49V09XPZydPnsTx48eRm5sLCwsLrdSnLpmZmWjbti0qKipqPZcYsn79+uHPP//U\ndTUMhnEdPXpAJBLh4MGDePz4Mc6dO4ezZ882+tcuEUGhUCArKwsdO3Z8qnpUVFQ81XJClpWVhb/9\n7W91njRbYtxikUjUIutpisbGnoiQnZ2NwMBALdWICUFT49zU1sjq57OsrCz4+fk9VQLWXOc3oX+W\nmR4gJih+fn6UmJiomn7zzTfpxRdfJCKiX375hXr37k329vbUrVs3Sk5OVs03YMAAWrx4MYWGhpKl\npSVNmjSJTE1NyczMjGxsbCgxMZHKy8tp7ty55OHhQR4eHjRv3jwqLy8nIqITJ06Qp6cnrVmzhtzd\n3SkqKopiY2Np7NixNGnSJJJIJNSlSxe6evUqrVy5klxdXcnHx4eOHTumqsOmTZuoY8eOJJFIqG3b\ntrRhwwbVe5Xlr127llxdXal169a0efNm1fulpaW0YMEC8vX1JTs7O+rbty+VlZU1uN3VpaWl0YAB\nA8je3p46depE3377LRERLV26lMzMzMjU1JRsbGxo06ZNNZZdtmwZTZo0STVd33rr29aCggJ64YUX\nyN7enhwdHalfv36kUCho0qRJJBaLydLSkmxsbOjjjz+uUYe6liUiOnfuHAUFBZFEIqHx48fT+PHj\n6b333iMios2bN1Pfvn3VyhKJRJSRkUFERAcPHqTu3buTra0teXt7U2xsrGq+mzdvkkgkori4OPLx\n8aEBAwYQEVFcXBx17NiRHBwcaMiQIZSVlVWjvk+ePCFra2sSiURkbW1NAQEBRETk6+tLx48fpyNH\njqjt9+7du9caN19fX9Vxv2zZskYddwMGDKBFixZRcHAw2dra0siRI6mwsFBt2+RyORERPXz4kKZO\nnUqtW7cmT09Peu+991Tvbd68mfr06UPz588ne3t78vf3p9OnT9OmTZvI29ubXF1d6auvvlLb9oUL\nF5KPjw+5ubnRzJkzVcdsfcf7hg0b1D6bI0aMICKiVatWkb+/P0kkEgoMDKR9+/YRkfKYtrCwIBMT\nE7KxsSEHBwciIoqOjlbFn4ho48aNFBAQQI6OjjRixAjKzc1VOxb++9//Urt27cje3p5mz56teu/a\ntWvUv39/srOzI2dnZxo/frzGca7r81ZZv5kzZ9KwYcPI2tpa7bxW6caNG9S/f3+SSCT03HPP0Zw5\nc1SfwcrYVVRUUHR0tNo+27Bhg9o+qTyev/vuO+rWrRvZ29tTnz596OLFi2rH2Jo1a6hLly5kYWFB\ncrm8wXPqkiVLKDQ0lCQSCQ0ePJju3btHRETe3t4kEonIxsaGbGxs6Ndff62xbWfOnKEePXqQra0t\nubm50YIFC1TvjR07ltzd3cnOzo769+9Pqampavvttddeo2HDhpGNjQ317duX8vLy6I033iB7e3vq\n0KED/fHHH2rbtWrVKgoMDCQHBweaMmUKPXnyhIiUx6GXl5dq3pycHBo9ejS5uLhQmzZtaN26dTXq\nXenQoUMUGBhIEomEPD096ZNPPqlR5s6dO1X7wMbGhszMzCgsLEx1zNT1+ajvPGdsOAkTGD8/Pzp+\n/DgREWVnZ1OnTp1o6dKldPv2bXJycqIjR44QEdEPP/xATk5OqpPCgAEDyNfXl9LS0kgul5NMJqOY\nmBhasmSJquwlS5ZQ7969qaCggAoKCqhPnz6q90+cOEGtWrWiRYsWkVQqpbKyMlq2bBlZWFjQsWPH\nqKKigiZPnky+vr60cuVKqqiooC+++ILatGmjKv/QoUN048YNIiL68ccfycrKis6dO6dW/rJly6ii\nooIOHz5MVlZW9PDhQyIimjVrFg0cOJByc3NVJ8fy8vI6t7ugoKDGvpNKpeTv70+rVq0imUxGSUlJ\nJJFIKD09nYiIYmNjKSoqqs59XzUJa2h/17atlSfGRYsW0cyZM6miooIqKiro1KlTavGt7cuoUl3L\nlpeXk4+PD3366adUUVFBe/fuJVNTU1X8GkrCkpOT6fLly0REdPHiRXJzc6P9+/cT0V9fdtHR0VRa\nWkplZWW0f/9+CggIoD///JPkcjl9+OGH1KdPnzrrXXVd1bezof1eff7GHncDBgwgT09PSk1NpZKS\nEhozZkyNL/LKRGvUqFE0c+ZMKi0tpfz8fAoODlYl0Js3b6ZWrVpRfHw8KRQKeu+998jT05PmzJlD\nUqmUjh07RhKJhEpKSoiIaN68eTRy5Eh68OABPX78mIYPH07vvPMOETV8vFf/bBIR7dmzh/Ly8oiI\naNeuXWRtbU137twhIqL4+Pga8a1aRmJiIjk7O9Mff/xB5eXl9Prrr1P//v3V4jN8+HAqKiqi7Oxs\ncnFxoaNHjxIR0YQJE2jlypVEpDzOTp8+XWecqsa5oc9bdHQ02dnZ0c8//0xEpEoMqgoJCaGFCxeS\nVCqln376iSQSiepYqR676vus+j45d+4cubq6UkpKCikUCvrqq6/Iz8+PpFIpESmTlaCgILp9+zY9\nefJEo3NqQEAAXbt2jcrKyigsLIwWLVpERESZmZlqdatNSEgIJSQkEBFRSUmJWqK2efNmKi4uJqlU\nSvPmzVP7cRIdHU3Ozs507tw5evLkCQ0aNIh8fX1p69atquNy4MCBqvl9fX2pS5cudPv2bSosLKTQ\n0FBVcl41YZLL5fTMM8/Q8uXLSSaT0Y0bN6ht27aq46A6d3d31fnn4cOHaufyqoldpUePHlHHjh1p\n48aNRFT/56O+c6Sx4SRMYHx9fcnGxobs7e3J19eXZs+eTWVlZbR69eoaX2RDhgxR/TIPCwujZcuW\nqb0fExOj9kvZ399fdcIhIjp69Cj5+fkRkfKDZWZmpmoZI1J+GQ4ePFg1/e2335KNjY3qF8ujR49I\nJBJRUVFRrdsyatQo+r//+z9V+ZaWlmonLVdXVzpz5gzJ5XKytLRU+9VaqaHtruqnn34id3d3tdcm\nTpyo+pVcvaWruqrvN2a91bd16dKlNHLkSLp+/XqN+RpKwupa9scffyQPDw+116om0Q0lYdXNnTuX\n5s+fT0R/fdndvHlT9f7QoUMpLi5ONS2Xy8nKyoqys7NrLa++JKyh/V7b/I057sLCwlQndyJl64yZ\nmRkpFAq1L/I7d+6Qubm56tc4EdH27dtVX2ibN2+mdu3aqd67ePEiiUQiys/PV73m5OREFy5cIIVC\nQdbW1mrb/PPPP6uSw/qOd6Ka
2014-05-05 05:44:10 +00:00
"text": [
2014-05-05 19:34:22 +00:00
"<matplotlib.figure.Figure at 0x10f9e1050>"
2014-05-05 05:44:10 +00:00
]
}
],
2014-05-05 19:34:22 +00:00
"prompt_number": 34
2014-05-05 05:44:10 +00:00
},
2014-05-04 22:57:02 +00:00
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name=\"cython_bonus\"></a>\n",
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bonus: How to use Cython without the IPython magic"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[[back to top](#sections)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"IPython's notebook is really great for explanatory analysis and documentation, but what if we want to compile our Python code via Cython without letting IPython's magic doing all the work? \n",
"These are the steps you would need."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 1. Creating a .pyx file containing the the desired code or function."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%file ccy_classic_lstsqr.pyx\n",
"\n",
"def ccy_classic_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" x_avg = sum(x)/len(x)\n",
" y_avg = sum(y)/len(y)\n",
" var_x = sum([(x_i - x_avg)**2 for x_i in x])\n",
" cov_xy = sum([(x_i - x_avg)*(y_i - y_avg) for x_i,y_i in zip(x,y)])\n",
" slope = cov_xy / var_x\n",
" y_interc = y_avg - slope*x_avg\n",
" return (slope, y_interc)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Writing ccy_classic_lstsqr.pyx\n"
]
}
],
"prompt_number": 11
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 2. Creating a simple setup file"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%file setup.py\n",
"\n",
"from distutils.core import setup\n",
"from distutils.extension import Extension\n",
"from Cython.Distutils import build_ext\n",
"\n",
"setup(\n",
" cmdclass = {'build_ext': build_ext},\n",
" ext_modules = [Extension(\"ccy_classic_lstsqr\", [\"ccy_classic_lstsqr.pyx\"])]\n",
")"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Writing setup.py\n"
]
}
],
"prompt_number": 12
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"####3. Building and Compiling"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"!python3 setup.py build_ext --inplace"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"running build_ext\r\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"cythoning ccy_classic_lstsqr.pyx to ccy_classic_lstsqr.c\r\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"building 'ccy_classic_lstsqr' extension\r\n",
"creating build\r\n",
"creating build/temp.macosx-10.6-intel-3.4\r\n",
"/usr/bin/clang -fno-strict-aliasing -Werror=declaration-after-statement -fno-common -dynamic -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -arch i386 -arch x86_64 -g -I/Library/Frameworks/Python.framework/Versions/3.4/include/python3.4m -c ccy_classic_lstsqr.c -o build/temp.macosx-10.6-intel-3.4/ccy_classic_lstsqr.o\r\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\u001b[1mccy_classic_lstsqr.c:2040:28: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyObject_AsString'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE char* __Pyx_PyObject_AsString(PyObject* o) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:2037:32: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function\r\n",
" '__Pyx_PyUnicode_FromString' [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE PyObject* __Pyx_PyUnicode_FromString(char* c_str) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:2104:26: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyObject_IsTrue'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE int __Pyx_PyObject_IsTrue(PyObject* x) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:2159:33: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyIndex_AsSsize_t'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE Py_ssize_t __Pyx_PyIndex_AsSsize_t(PyObject* b) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:2188:33: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyInt_FromSize_t'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE PyObject * __Pyx_PyInt_FromSize_t(size_t ival) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:1584:32: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyInt_From_long'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE PyObject* __Pyx_PyInt_From_long(long value) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:1631:27: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1mfunction '__Pyx_PyInt_As_long' is not\r\n",
" needed and will not be emitted [-Wunneeded-internal-declaration]\u001b[0m\r\n",
"static CYTHON_INLINE long __Pyx_PyInt_As_long(PyObject *x) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:1731:26: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1mfunction '__Pyx_PyInt_As_int' is not\r\n",
" needed and will not be emitted [-Wunneeded-internal-declaration]\u001b[0m\r\n",
"static CYTHON_INLINE int __Pyx_PyInt_As_int(PyObject *x) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"8 warnings generated.\r\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\u001b[1mccy_classic_lstsqr.c:2040:28: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyObject_AsString'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE char* __Pyx_PyObject_AsString(PyObject* o) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:2037:32: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function\r\n",
" '__Pyx_PyUnicode_FromString' [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE PyObject* __Pyx_PyUnicode_FromString(char* c_str) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:2104:26: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyObject_IsTrue'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE int __Pyx_PyObject_IsTrue(PyObject* x) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:2159:33: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyIndex_AsSsize_t'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE Py_ssize_t __Pyx_PyIndex_AsSsize_t(PyObject* b) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:2188:33: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyInt_FromSize_t'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE PyObject * __Pyx_PyInt_FromSize_t(size_t ival) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:1584:32: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1munused function '__Pyx_PyInt_From_long'\r\n",
" [-Wunused-function]\u001b[0m\r\n",
"static CYTHON_INLINE PyObject* __Pyx_PyInt_From_long(long value) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:1631:27: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1mfunction '__Pyx_PyInt_As_long' is not\r\n",
" needed and will not be emitted [-Wunneeded-internal-declaration]\u001b[0m\r\n",
"static CYTHON_INLINE long __Pyx_PyInt_As_long(PyObject *x) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m\u001b[1mccy_classic_lstsqr.c:1731:26: \u001b[0m\u001b[0;1;35mwarning: \u001b[0m\u001b[1mfunction '__Pyx_PyInt_As_int' is not\r\n",
" needed and will not be emitted [-Wunneeded-internal-declaration]\u001b[0m\r\n",
"static CYTHON_INLINE int __Pyx_PyInt_As_int(PyObject *x) {\r\n",
"\u001b[0;1;32m ^\r\n",
"\u001b[0m"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"8"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
" warnings generated.\r\n",
"/usr/bin/clang -bundle -undefined dynamic_lookup -arch i386 -arch x86_64 -g build/temp.macosx-10.6-intel-3.4/ccy_classic_lstsqr.o -o /Users/sebastian/Github/python_reference/benchmarks/ccy_classic_lstsqr.so\r\n"
]
}
],
"prompt_number": 13
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 4. Importing and running the code"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import ccy_classic_lstsqr\n",
"\n",
"%timeit classic_lstsqr(x, y)\n",
"%timeit cy_classic_lstsqr(x, y)\n",
"%timeit ccy_classic_lstsqr.ccy_classic_lstsqr(x, y)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"100 loops, best of 3: 2.9 ms per loop\n",
"1000 loops, best of 3: 212 \u00b5s per loop"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"1000 loops, best of 3: 207 \u00b5s per loop"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n"
]
}
],
"prompt_number": 20
2014-05-05 19:34:22 +00:00
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name=\"numba\"></a>\n",
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Appendix I: Cython vs. Numba"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[[back to top](#sections)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Like we did with Cython before, we will use the minimalist approach to Numba and see how they compare against each other. \n",
"\n",
"Numba is using the [LLVM compiler infrastructure](http://llvm.org) for compiling Python code to machine code. Its strength is to work with NumPy arrays to speed-up code. If you want to read more about Numba, see the original [website and documentation](http://numba.pydata.org/numba-doc/0.13/index.html)"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" x_avg = sum(x)/len(x)\n",
" y_avg = sum(y)/len(y)\n",
" var_x = sum([(x_i - x_avg)**2 for x_i in x])\n",
" cov_xy = sum([(x_i - x_avg)*(y_i - y_avg) for x_i,y_i in zip(x,y)])\n",
" slope = cov_xy / var_x\n",
" y_interc = y_avg - slope*x_avg\n",
" return (slope, y_interc)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%load_ext cythonmagic"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%cython\n",
"\n",
"def cy_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" x_avg = sum(x)/len(x)\n",
" y_avg = sum(y)/len(y)\n",
" var_x = sum([(x_i - x_avg)**2 for x_i in x])\n",
" cov_xy = sum([(x_i - x_avg)*(y_i - y_avg) for x_i,y_i in zip(x,y)])\n",
" slope = cov_xy / var_x\n",
" y_interc = y_avg - slope*x_avg\n",
" return (slope, y_interc)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from numba import jit\n",
"\n",
"@jit\n",
"def nmb_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" x_avg = sum(x)/len(x)\n",
" y_avg = sum(y)/len(y)\n",
" var_x = sum([(x_i - x_avg)**2 for x_i in x])\n",
" cov_xy = sum([(x_i - x_avg)*(y_i - y_avg) for x_i,y_i in zip(x,y)])\n",
" slope = cov_xy / var_x\n",
" y_interc = y_avg - slope*x_avg\n",
" return (slope, y_interc)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"ename": "ImportError",
"evalue": "No module named 'llvm'",
"output_type": "pyerr",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mImportError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-4-505a644bdc72>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mnumba\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mjit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mjit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mnmb_lstsqr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.3/lib/python3.3/site-packages/numba/__init__.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 3\u001b[0m \"\"\"\n\u001b[1;32m 4\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0m__future__\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mprint_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdivision\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mabsolute_import\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0;34m.\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mtesting\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdecorators\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 6\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0;34m.\u001b[0m\u001b[0m_version\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mget_versions\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0;31m# Re-export typeof\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.3/lib/python3.3/site-packages/numba/decorators.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mwarnings\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mnumba\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0msigutils\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mnumba\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtargets\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mregistry\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;31m# -----------------------------------------------------------------------------\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.3/lib/python3.3/site-packages/numba/targets/registry.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0m__future__\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mprint_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdivision\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mabsolute_import\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mnumba\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mutils\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtyping\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mnumba\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtargets\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mcpu\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mnumba\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtargets\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdescriptors\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mTargetDescriptor\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mnumba\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mdispatcher\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.3/lib/python3.3/site-packages/numba/targets/cpu.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mllvm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcore\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mlc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 5\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mllvm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpasses\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mlp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mllvm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mee\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mle\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mImportError\u001b[0m: No module named 'llvm'"
]
}
],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ... this section is still in progress"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name=\"type_declarations\"></a>\n",
"<br>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Appendix II: Cython with and without type declarations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[[back to top](#sections)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the sections above, we have been using the simplest approach to Cython without using static type declarations and thereby neglecting one of its major strengths. \n",
"Let us now see how we can further improve the Cython implementation of our \"classic least squares approach\" by adding those type declarations."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, our \"simple\" approach using Cython from the previous section:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%cython\n",
"def cy_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" x_avg = sum(x)/len(x)\n",
" y_avg = sum(y)/len(y)\n",
" var_x = sum([(x_i - x_avg)**2 for x_i in x])\n",
" cov_xy = sum([(x_i - x_avg)*(y_i - y_avg) for x_i,y_i in zip(x,y)])\n",
" slope = cov_xy / var_x\n",
" y_interc = y_avg - slope*x_avg\n",
" return (slope, y_interc)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 54
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And now, the same code with static type declarations:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%%cython\n",
"def static_type_lstsqr(x, y):\n",
" \"\"\" Computes the least-squares solution to a linear matrix equation. \"\"\"\n",
" cdef double x_avg, y_avg, var_x, cov_xy, slope, y_interc, x_i, y_i\n",
" x_avg = sum(x)/len(x)\n",
" y_avg = sum(y)/len(y)\n",
" var_x = sum([(x_i - x_avg)**2 for x_i in x])\n",
" cov_xy = sum([(x_i - x_avg)*(y_i - y_avg) for x_i,y_i in zip(x,y)])\n",
" slope = cov_xy / var_x\n",
" y_interc = y_avg - slope*x_avg\n",
" return (slope, y_interc)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 55
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>\n",
"Now, let us see how the two functions (with and without static type declarations) compare against each other for different sample sizes."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import timeit\n",
"import random\n",
"random.seed(12345)\n",
"\n",
"funcs = ['cy_lstsqr', 'static_type_lstsqr'] \n",
"labels = ['simple Cython', 'Cython w. type declarations']\n",
"orders_n = [10**n for n in range(1, 7)]\n",
"times_n = {f:[] for f in funcs}\n",
"\n",
"for n in orders_n:\n",
" x = [x_i*random.randrange(8,12)/10 for x_i in range(n)]\n",
" y = [y_i*random.randrange(10,14)/10 for y_i in range(n)]\n",
" for f in funcs:\n",
" times_n[f].append(timeit.Timer('%s(x,y)' %f, \n",
" 'from __main__ import %s, x, y' %f).timeit(1000))"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 58
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#%pylab inline\n",
"#import matplotlib.pyplot as plt\n",
"\n",
"plt.figure(figsize=(10,8))\n",
"\n",
"for f in times_n.keys():\n",
" plt.plot(orders_n, times_n[f], alpha=0.5, label=f, marker='o', lw=2)\n",
"\n",
"plt.xlabel('sample size n')\n",
"plt.ylabel('time in ms')\n",
"plt.xlim([0,max(orders_n) + max(orders_n) * 0.1])\n",
"plt.legend(loc=2)\n",
"plt.grid()\n",
"\n",
"plt.title('Performance of a simple least square fit implementation')\n",
"plt.show()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAH4CAYAAAACdDpdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYVGf6N/Dv0KQN0osComhUXBVs2BJRsUaxYzCiJkZT\nNm3NrtFsVEyylrzJpq1J1KhYY4sNW4gFW1xL1JiIgqKAAqIoFqQMwzzvH+fnLChNYZgzc76f6+K6\n5szMOdwz9wxz8zz3PEclhBAgIiIiojplYewAiIiIiJSIRRgRERGREbAIIyIiIjICFmFERERERsAi\njIiIiMgIWIQRERERGQGLMDI52dnZeO655+Dk5IR//OMfxg7H6AoKCjB48GA4Oztj9OjRBvs9AwcO\nxMqVKw1ybAsLC1y+fPmJ94uNjcWzzz5rgIjM35EjR9CsWTM4OTlh69atGDhwIFasWFHt/eX4epC7\nuXPnYtKkScYOg2TEytgBkDIEBATgxo0bsLS0hIODAwYMGID//Oc/cHBweOJjLVq0CJ6enrh3754B\nIjU9GzduxI0bN3D79m1YWBju/6qdO3ca7NhyFhYWhujoaEycONHYodSqmTNn4u2338Zbb70FABgy\nZIj+ttjYWCxZsgSHDh2qcH+lvR5iYmKQkpJS7cIzISEB0dHRuHr1qv666dOnGyo8MlEcCaM6oVKp\nsH37dty/fx+nTp3CyZMn8cknnzzRMYQQ0Ol0SEtLQ8uWLZ8qDq1W+1T7yVlaWhqeeeYZgxZgSqZS\nqYwdQpWe5nWdnp6OoKAgA0RDRNUmiOpAQECA2Lt3r37773//uxg0aJAQQoijR4+KLl26CGdnZ9G2\nbVuRkJCgv1+PHj3EP//5T9GtWzdhZ2cnxo4dK6ytrYWNjY1wdHQUe/fuFUVFReKdd94RDRo0EA0a\nNBDvvvuuKCoqEkIIsX//ftGwYUMxf/584e3tLaKjo0VMTIwYOXKkGDt2rFCr1aJ169YiOTlZzJkz\nR3h6egp/f38RHx+vj2Hp0qWiZcuWQq1WiyZNmoiFCxfqb3t4/M8//1x4enoKHx8fsWzZMv3t+fn5\nYsqUKaJRo0aifv36onv37qKgoKDKx/2oxMRE0aNHD+Hs7CxatWoltm3bJoQQYubMmcLGxkZYW1sL\nR0dHsXTp0sf2PXbsmOjcubNwdnYWPj4+4s033xQajabc31NQUCBefPFF4ebmJpydnUXHjh3FjRs3\n9Ln44YcfhBBCLFu2THTt2lX87W9/E87OziIwMFAcOXJELF26VPj5+QlPT0+xfPly/XHHjx8vXn31\nVdGnTx+hVqtFjx49RFpamv52lUolUlJShBBCFBYWivfee0/4+/sLLy8v8dprr+mfs0ctW7ZMdO/e\nXb99/vx5ER4eLlxdXUXz5s3F+vXr9bdt375dBAcHCycnJ+Hn5ydiYmIqfdzZ2dnigw8+EJaWlsLW\n1lY4OjqKt956q1rPWXZ2thBCiMuXL4vnnntOqNVq0adPH/HXv/5VjB07VgghvXZ8fX3LHKtRo0b6\n90lVeVOpVGLBggWiadOmokmTJkIIIeLi4kTbtm2Fs7Oz6Nq1qzh79my5z1uTJk2EhYWFsLOzE2q1\nWhQVFenze/78eVGvXj1haWkpHB0dhYuLS7nHkMPr4eH779NPPxUeHh7Cx8dHbN68WezYsUM0a9ZM\nuLq6irlz5+qPq9PpxNy5c0VgYKBwc3MTkZGR4vbt20IIIa5cuSJUKpVYvny58Pf3F+7u7uJf//qX\nEEKIXbt2lXmfBQcHCyEq/tuQl5cnbG1thYWFhXB0dBRqtVpkZmaKWbNm6fMvhBBbt24VQUFBwtnZ\nWYSFhYnz58+XeS189tlnok2bNqJ+/fpi9OjRorCwsNxckOliEUZ1IiAgQOzZs0cIIUR6erpo1aqV\nmDlzprh27Zpwc3MTu3btEkII8csvvwg3NzeRk5MjhJD+0Ddq1EgkJiaKkpISUVxcLCZMmCBmzJih\nP/aMGTNEly5dxM2bN8XNmzdF165d9bfv379fWFlZiWnTpgmNRiMKCgrErFmzhK2trYiPjxdarVaM\nGzdONGrUSMyZM0dotVqxePFi0bhxY/3xd+zYIS5fviyEEOLAgQPC3t5enDp1qszxZ82aJbRardi5\nc6ewt7cXd+7cEUII8cYbb4iePXuKzMxMUVJSIo4ePSqKiooqfNw3b9587LnTaDQiMDBQzJ07VxQX\nF4t9+/YJtVotkpKShBBCxMTEiOjo6Aqf+99++00cO3ZMlJSUiNTUVNGyZUvx5Zdflnvf77//Xgwe\nPFgUFBQInU4nTp06Je7duyeEECIsLEwsWbJECCF96FpZWYnY2Fih0+nEhx9+KBo2bKgvFOLj44Va\nrRYPHjwQQkgfumq1Whw6dEhfNJcunkp/6L777rtiyJAhIjc3V9y/f18MHjxYTJ8+vdx4SxdheXl5\nwtfXV8TGxoqSkhJx+vRp4e7uLhITE4UQQiQkJIg///xTCCHE2bNnhZeXl9iyZcsTPe4nfc46d+4s\n3nvvPaHRaMTBgweFWq3W56q8Iqz0PytV5U2lUom+ffuK3NxcUVhYKE6dOiU8PT3F8ePHhU6nE8uX\nLxcBAQH6f0ge9eg/RqUfZ2xsbJn8lEcOr4eH77+PP/5Y/951c3MTY8aMEXl5eeLcuXPCzs5OpKam\nCiGE+PLLL0WXLl1ERkaG0Gg04tVXXxVRUVFCiP8VYZMnTxaFhYXi999/F/Xq1RMXLlwQQpT/Pqvs\nb0NCQsJj+Y2JidEXYUlJScLBwUHs2bNHaLVa8emnn4qmTZuK4uJifX5CQ0NFVlaWuH37tmjZsqX4\n/vvvK80JmR4WYVQnGjVqJBwdHYWzs7No1KiR+Otf/yoKCgrEvHnzHvvD1q9fP/1/zWFhYWLWrFll\nbp8wYYL48MMP9duBgYH6YkYIIX7++WcREBAghJD+SNvY2JT5IJo1a5bo27evfnvbtm3C0dFR6HQ6\nIYQQ9+7dEyqVSty9e7fcxzJ06FDx1Vdf6Y9vZ2cnSkpK9Ld7enrqPzzt7OzKHY2o6nGXdvDgQeHt\n7V3muqioKP1IzqP/XVfliy++EMOGDSv3tqVLl1Y4gvLoh26zZs30t509e1aoVCr9qJkQQri5uYnf\nf/9dCCF96D78sBNCKpgsLS3FtWvXhBD/+9DV6XTCwcFB/wEshBC//vprmaK4tNJF2Nq1a8Wzzz5b\n5vbJkyeL2bNnl7vvO++8I/72t79V63E/HPEpT0X7pqWlCSsrK5Gfn6+/bsyYMdUuwh71aN5UKpXY\nv3+/fvu1114r88+JEEI0b95cHDhwoNzjVVaEPTrCWB45vB4evv8efe8eP35cf//27duLrVu3CiGE\naNGiRZnHnJmZKaytrUVJSYm+CMvIyNDf3qlTJ7Fu3TohRPXeZ4/+bXg0v6WP8dFHH4nRo0frb9Pp\ndKJhw4b6fAUEBIjVq1frb586dap47bXXKv39ZHrYmE91QqVSYevWrejVq1eZ69PS0rBhwwbExcXp\nr9NqtWXu5+fnV+mxMzMz0ahRI/22v78/MjMz9dseHh6wsbEps4+np6f+sp2dHdzd3fW9P3Z2dgCA\nvLw8ODk5YdeuXZg9ezYuXrwInU6H/Px8tGnTRr+/m5tbmX4se3t75OXlIScnB4WFhQgMDHws5uo8\n7tKP79HnoFGjRsjIyKj0eXkoOTkZU6ZMwW+//Yb8/HxotVp06NCh3Ps+bCR+4YUXcOfOHYwdOxb/\n+te/YGX1+J8KLy8v/eWHz5mHh0eZ6/Ly8gBI+ff19dXf5uDgAFdXV2RmZqJhw4b662/evIn8/Hy0\nb99ef534v17AqqSlpeHYsWNwcXHRX6fVajFu3DgAwLFjxzBt2jScO3cOGo0GRUVFiIyMrNbjrqwv\nrKJ9MzMz4eLion9uAClvpRu1K1OdvJV+XaSlpWHFihX45ptv9NcVFxcjKyurWr+vpoz1enBzc3vs\nvftoLA9/b1paGoYNG1bm/Wpl
"text": [
"<matplotlib.figure.Figure at 0x10eea0790>"
]
}
],
"prompt_number": 59
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<br>\n",
"<br>\n",
"The improvement is pretty significant when static type declarations are used. One more experiment to see by how much we could improve our \"classic least squares\" code via Cython compared to the initial Python implementation.\n",
"<br>\n",
"<br>"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import random\n",
"random.seed(12345)\n",
"\n",
"x = [x_i*random.randrange(8,12)/10 for x_i in range(500)]\n",
"y = [y_i*random.randrange(8,12)/10 for y_i in range(100,600)]"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 60
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import timeit\n",
"\n",
"funcs = ['classic_lstsqr', 'cy_lstsqr', 'static_type_lstsqr', \n",
" 'lin_lstsqr_mat', 'numpy_lstsqr', 'scipy_lstsqr']\n",
"labels = ['classic approach','classic approach (cython)', \n",
" 'classic approach (cython + type decl.)',\n",
" 'matrix approach', 'numpy function', 'scipy function']\n",
"\n",
"times = [timeit.Timer('%s(x,y)' %f, \n",
" 'from __main__ import %s, x, y' %f).timeit(1000)\n",
" for f in funcs]\n",
"\n",
"times_rel = [times[0]/times[i+1] for i in range(len(times[1:]))]"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 61
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#%pylab inline\n",
"#import matplotlib.pyplot as plt\n",
"\n",
"x_pos = np.arange(len(funcs[1:]))\n",
"plt.bar(x_pos, times_rel, align='center', alpha=0.5, color=\"green\")\n",
"plt.xticks(x_pos, labels[1:], rotation=90)\n",
"plt.ylabel('relative performance gain')\n",
"plt.title('Performance gain compared to the classic least square implementation')\n",
"plt.grid()\n",
"plt.show()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAbwAAAG/CAYAAAAjJGJxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8Dff3P/DXJGLLLhFJSHIRSgiJKLEmqb3WWoLakqBo\naX3ThSoV9aH0Ux+EVmk1UdUWpYSiisYeEhVbKrFdSyQRImSP3Lx/f+SXaa5ESNw7c+d9z/PxyOOR\nu82cc+/cOXfmzLxHYIwxEEIIIZwzkTsAQgghRApU8AghhBgFKniEEEKMAhU8QgghRoEKHiGEEKNA\nBY8QQohRkLTgpaWloXv37rCyssKHH34o5awV6fXXX8fGjRvlDsNgqVQqHDx4ULL5mZiY4Pr163qd\nR1BQEObNm6e36VtaWkKtVlfpNdHR0XBxcdFPQJxr3bo1jhw5ovPpqtVqmJiYoLi4WOfTltu0adPw\nn//8Ry/TrvG8J6hUKty7dw+mpqYwNzdHv379sHr1apibm1d5ZuvWrYODgwMeP35crWCNzZ49e+QO\nwaAJggBBECp8LCgoCC4uLli4cGG1pu3v749x48Zh4sSJLxNilVWWky5kZWXpbdov42U/L0N18eJF\nuUOQVFU/x8jISKxfvx5Hjx4V71uzZo2+wnv+Fp4gCNi9ezeysrLw999/Iy4ursrVlzGG4uJi3Lx5\nEy1btqxWoEVFRdV6HVEGQ/t89Vl0nofGgjBMhraMkmpgz6FSqdjBgwfF2x988AEbMGAAY4yxkydP\nsk6dOjEbGxvWtm1bFh0dLT7Pz8+PffLJJ6xLly6sTp06bOzYsczMzIzVrFmTWVhYsIMHD7KCggL2\n3nvvMWdnZ+bs7MxmzpzJCgoKGGOM/fXXX6xhw4Zs6dKlzNHRkY0bN46FhYWx4cOHs7FjxzJLS0vm\n6enJkpKS2OLFi5mDgwNzdXVl+/fvF2P4/vvvWcuWLZmlpSVr0qQJW7t2rfhY6fSXLVvGHBwcmJOT\nE4uIiBAfz83NZaGhoczNzY1ZW1uzrl27sry8vOfm/bQzZ84wLy8vZmlpyUaMGMECAwPZ3LlzGWOM\nZWRksP79+7P69eszW1tbNmDAAHbnzh2t9/C7775jjDEWERHBunTpwj744ANma2vLGjduzPbu3fvM\n+d66dYu98cYbrH79+szOzo5Nnz6dMcaYRqNhCxcuZG5ubszBwYGNHz+ePXr0iDHG2I0bN5ggCCwi\nIoK5uLiwevXqsTVr1rDTp08zT09PZmNjI06nNKbOnTuz6dOnM2tra9aiRQutZeVF3v/Sz3f8+PGs\nuLiYff7556xp06bMzs6OBQYGsoyMDPE1P/zwA3N1dWV2dnZs0aJF5ZbNUmvXrtVa1gYNGsQYYywh\nIYH5+fkxGxsb1qpVKxYVFVXhezdnzhxmamrKateuzSwsLNiMGTMYY4wJgsC++eYb1qxZM2ZjY8Pe\neecdrdetX7+etWzZktna2rI+ffqwmzdvPvPzOXr0qLgMubi4sA0bNjDGGAsKCnrh5SMiIoI1adKE\nWVpassaNG7NNmzYxxhi7cuUK6969O7O2tmb29vZs5MiR4msEQWDXrl1jjFW+jJf1119/sUaNGom3\nk5OT2dChQ1n9+vVZ48aNWXh4uPjYqVOnmK+vL7OxsWFOTk5s+vTprLCwUHx85syZzMHBgVlZWTFP\nT0928eLFZ35eT6votYwxdv/+fTZw4EBmZWXFOnTowObOncu6du3KGPt3mdZoNOJ0yn6vrl69ygIC\nApidnR2zt7dnY8aMYZmZmeJz3dzc2NKlS5mnpyerXbs202g0Vfr+u7m5icvo/Pnzq7T+8vPzY7Nn\nz2YdOnRgVlZWbPDgweL34em8MjMzWUhICHNycmINGzZkc+fOFR8r/Z7+3//9H7OxsWFNmzZlx48f\nZ99//z1zcXFhDg4O4vLHGGP5+fns/fffZ66urqxBgwZs6tSp4nJR2XrzWZ9j6Xfa0tKSeXh4sN9+\n+40xVvJ9rF27NjM1NWUWFhbM1taWMcbYhAkTxO8AY4ytW7eOubu7s3r16rFBgwaxu3fvio897zv5\ntBcqeAcOHGCMlaxEW7VqxT799FN2584dZmdnJ650//zzT2ZnZ8fu378vflhubm4sISGBaTQa9uTJ\nExYUFMTmzZsnTnvevHmsU6dOLD09naWnp7POnTuLj//111+sRo0abPbs2aywsJDl5eWx+fPns9q1\na7P9+/ezoqIiNn78eObm5sYWL17MioqK2LfffssaN24sTv/3339n169fZ4wxdvjwYVa3bl32999/\na01//vz5rKioiO3Zs4fVrVtXXNjffvttFhAQwO7evSsu5AUFBc/MOz09vdx7V1BQwFxdXVl4eDgr\nKipi27dvZzVr1hRzfPDgAdu+fTvLy8tjWVlZbMSIEWzIkCHi6/39/dn69esZYyULrZmZGfvuu+9Y\ncXExW7NmDXN2dq7wMysqKmJt2rRhoaGhLDc3l+Xn57Pjx48zxkpWyu7u7uzGjRssOzubDR06lI0b\nN44x9u+XaNq0aaygoIDt37+f1axZkw0ZMoSlp6ez5ORk5uDgwA4fPizGVKNGDbZixQpWVFTENm/e\nzKytrcUv5Yu8/2U/3xUrVrBOnTqx5ORkVlhYyKZMmcJGjx7NGGPs0qVLzMLCgh09epQVFBSw0NBQ\nVqNGjQoLHmOs3LJWWFjImjZtyj7//HP25MkTdujQIWZpackSExMrfH3Z976UIAhs4MCB7NGjR+zW\nrVusfv36bN++fYwxxnbs2MHc3d3Z5cuXmUajYf/5z39Y586dK5y2Wq1mlpaW7JdffmFFRUXswYMH\nLD4+Xoy79Mte2fKRnZ3NrKysWFJSEmOMsdTUVHbp0iXGGGOjRo1iixcvZoyVLIOln31pDqUF71nL\n+NPKFjyNRsPatWvHFi5cyJ48ecKuX7/OmjRpwv744w/GWMkPvFOnTjGNRsPUajVr2bIlW7FiBWOM\nsX379jEfHx/xB9bly5dZSkpKhZ/X0yp77ciRI9nIkSNZbm4uu3jxImvYsCHr1q0bY6ziglf2s716\n9So7cOAAKywsZOnp6ax79+5s5syZ4nPd3NyYt7c3u3PnDsvPz6/S958x7Q2Gqq6//Pz8WMOGDdml\nS5dYTk4OGzZsGBs7dmyFeQ0ZMoRNnTqV5ebmsnv37rEOHTqIPzBLv6eRkZGsuLiYzZ07lzVs2FD8\nMbJ//35maWnJcnJyGGMlPywGDx7MHj58yLKystjAgQPZxx9/LC4Lla03K/oct27dKn5WmzdvZubm\n5iw1NZUxxlhkZKT446RU2WkcPHiQ2dvbs7Nnz7KCggI2Y8YM1r17d/G5lX0nK/Lcgufm5sYsLCyY\njY0Nc3NzY++88w7Ly8tjS5YsEVeUpfr06SP+UvD392fz588vl0jZyt20aVOtrZQ//viDqVQqxljJ\nG1uzZk2tL+D8+fNZ7969xdtRUVHMwsKCFRcXM8YYe/z4MRMEQfxSPG3IkCFs5cqV4vTr1Kmj9UVw\ncHAQv6x16tRh58+fLzeN5+Vd1uHDh1nDhg217uvateszv9hnz54Vf+UwVr7gubu7i4/l5OQwQRBY\nWlpauemcOHGC1a9fXyu3Uq+99hpbs2aNeDsxMZGZmZkxjUYjfonK/oKys7NjW7ZsEW8PGzZMXIFF\nRESUK7odOnRgGzdurDC/p9//pz/fli1bahWwu3fvMjMzM1ZUVMQWLFggFr/S/GvWrFlpwSu7rB05\ncoQ5OjpqPWf06NEsLCyswtf7+/uLWwGlBEHQKh6BgYFs6dKljDHG+vbtq1UgNRoNq1u3Lrt161a5\naS9evJgNHTr0heIuq+zykZ2dzWxsbNi2bdtYbm6u1vPGjx/P3nrrLa2twbI5XLt2rdJl/GllC15M\nTAxzdXUtl09wcHCFr12+fDl74403GGMlK6/mzZuzmJiYcstmZXkzxtihQ4cqfG1RUREzMzPT+uEy\nZ86cSrfwKvoxU+q3335j3t7e
"text": [
"<matplotlib.figure.Figure at 0x106418090>"
]
}
],
"prompt_number": 63
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a pretty significant performance gain. The \"Cython + type declarations\" approach sped up our initial Python code 25 times."
]
2014-05-04 22:57:02 +00:00
}
],
"metadata": {}
}
]
}
